/*
 *  process.cpp
 *  EPPI-0.0
 *
 *  Created by Juan Carlos Garcia Orden on 15/11/07.
 *	Modified by Markus Lutz
  *	Modified by Sergio Conde on 17/06/2010
 *
 */

#include "process.h"

Process::Process() 
:	num(0),
	lrank(rankk),
	dimension(3),
	mMecSystems(),
	rootSys(0),
	mSysProc(),
	vPTProc(),
	relatedProc(),
	vInterfaces(),
	localSharedPoints(),
//	lastAnalysis(),
//	currentAnalysis(),
	numDofs(0),
	q(),
	dq(),
	ddq(),
	qConv(),
	dqConv(),
	res(),
	Deltaq(),
	isVelNeededLast(false),
	mMaterials(),
	outFileName(),
	ou(),
	localGlobI(),
	globalLocalI(),
	localGlobB(),
	globalLocalB(),
	numDofsI(0),
	numDofsB(0),
	Kint(),
	vDofsOut(),
	vStressOut(),
	outDofsFile(),
	outVelDofsFile(),
	outGraphFile(),
	outTempFile(),
	outIterFile(),
	outStressFile(),
	fDofs(),
	fVelDofs(),
	fGraph(),
	fTemp(),
	fIter(),
	fStress(),
	vecType(0),
	matType(0),
	procSS(),
	grav(0.0),
	gravDir(),
	eneResFi(1.0e9),
	eneResFq(1.0e9),
	eneResFt(1.0e9),
	eneResInter(1.0e9),
	eneResMax(),
	eneRes(),
	resFq(),
	resFt(),
	resInter(),
	sparseIntTot(),
	M(),
	dynFlag(0),
	Ekin(0.0),
	Epot(0.0),
	Econ(0.0),
	Eext(0.0),
	Eamb(0.0),
	Dint(0.0),
	Entropy(0.0),
	EntropyAmb(0.0),
	AM(),
	vCon(),
	numDofsInertia(),
	numNoInertia(),
	thetab(300.0),
	vTempOut(),
	projFlag(0),
	iaInertia(),
	jaInertia(),
	projVel(),
	projPenaltyMax(1.0e8),
	Gamma(),
	GammaFac(),
	auxGamma(),
	symmFlag(true),
	timeInfoFileName(),
	ouTime()
{}

// ***

Process::Process(int n, int dim, int vecT, int matT, int& mode, const char *simTitle)
:	num(n),
	lrank(rankk),
	dimension(dim),
	mMecSystems(),
	rootSys(0),
	mSysProc(),
	vPTProc(),
	relatedProc(),
	vInterfaces(),
	localSharedPoints(),
//	lastAnalysis(),
//	currentAnalysis(),
	numDofs(0),
	q(),
	dq(),
	ddq(),
	qConv(),
	dqConv(),
	res(),
	Deltaq(),
	isVelNeededLast(false),
	mMaterials(),
	outFileName(),
	ou(),
	localGlobI(),
	globalLocalI(),
	localGlobB(),
	globalLocalB(),
	numDofsI(0),
	numDofsB(0),
	Kint(),
	vDofsOut(),
	vStressOut(),
	outDofsFile(),
	outVelDofsFile(),
	outGraphFile(),
	outTempFile(),
	outIterFile(),
	outStressFile(),
	fDofs(),
	fVelDofs(),
	fGraph(),
	fTemp(),
	fIter(),
	fStress(),
	vecType(vecT),
	matType(matT),
	procSS(),
	grav(0.0),
	gravDir(),
	eneResFi(1.0e9),
	eneResFq(1.0e9),
	eneResFt(1.0e9),
	eneResInter(1.0e9),
	eneResMax(),
	eneRes(),
	resFq(),
	resFt(),
	resInter(),
	sparseIntTot(),
	M(),
	dynFlag(0),
	Ekin(0.0),
	Epot(0.0),
	Econ(0.0),
	Eext(0.0),
	Eamb(0.0),
	Dint(0.0),
	Entropy(0.0),
	EntropyAmb(0.0),
	AM(),
	vCon(),
	numDofsInertia(),
	numNoInertia(),
	thetab(300.0),
	vTempOut(),
	projFlag(0),
	iaInertia(),
	jaInertia(),
	projVel(),
	projPenaltyMax(1.0e8),
	Gamma(),
	GammaFac(),
	auxGamma(),
	symmFlag(true),
	timeInfoFileName(),
	ouTime()
{
	AM.resize(dimension);
	vCon.resize(3);
	
	strcpy(title,simTitle);
	
	if(mode==0)	// General info
	{
		stringstream stnum;
		string snum; 
		stnum << num;
		stnum >> snum;
	
		strcpy(outFileName,"out-");
		strcat(outFileName,snum.c_str());
	
		stnum.clear();
	
		stnum << rankk;
		stnum >> snum;
		
		strcat(outFileName,".dat");

		if(num==rankk)
		{
			ou.open(outFileName);
			ou.close();
		}
	}
	if(mode==2 || mode==0)	// Time info
	{
		tini = MPI::Wtime();	// Starting time
		
		stringstream stnum;
		string snum; 
		stnum << num;
		stnum >> snum;
		
		strcpy(timeInfoFileName,"timeInfo-");
		strcat(timeInfoFileName,simTitle);
		strcat(timeInfoFileName,"-");
		strcat(timeInfoFileName,snum.c_str());
		
		stnum.clear();
		
		stnum << rankk;
		stnum >> snum;
		
		strcat(timeInfoFileName,".dat");
		
		if(num==rankk)
		{
			ouTime.open(timeInfoFileName);
			ouTime.close();
		}
	}
}

// ***

Process::~Process()
{
	if(ou.is_open())	ou.close();
	if(ouTime.is_open())	ouTime.close();
}

// ***

int Process::readSystem(string nam, ifstream &in, int dimension)
{
	mMecSystems[nam]= new MecSystem(nam, in,dimension,this,ou);
	if (mMecSystems[nam]->isWellDefined) return(1) ; else return (0);	
}

// ***

void Process::setSysProc(map<string, int, less<string> > msp, vector<Process*> vptp)
{
	mSysProc = msp;
	vPTProc = vptp;
}

// ***

Process* Process::getPTProc(string sysName)
{
	return( vPTProc[ mSysProc[sysName] ] );
}

// ***

Process* Process::getPTProc(int num)
{
	return( vPTProc[num] );
}

// ***

MecSystem* Process::getPTSys(Process* pp, string sysName)
{
	return( ( pp->mMecSystems[sysName] ) ); 
}

// ***

MecSystem* Process::getPTSystem(string sysName)
{
	return( getPTSys( getPTProc(sysName), sysName) ); 
}

// ***

void Process::buildSysTree()
{
	map<string, MecSystem*, less<string> >::iterator i;
	vector<MecSystem*> vsaux;	/** Auxiliar vector of orphan systems */

	if (num==rankk) cout << "	Building system tree in process " << num << " ...\n";
	/** Set parenthood relations among systems of this process */	
	for(i=mMecSystems.begin(); i!=mMecSystems.end(); i++)	((*i).second)->setParenthood();

	/** Set root system of the process */
	/** Begin loop on systems, checking for orphans */
	for(i=mMecSystems.begin(); i!=mMecSystems.end(); i++)
	{
		if( ((*i).second)->parentSysName=="" )
		{
			vsaux.push_back( (i->second) );
		}
	}
	/** Loop on systems completed */
	/** Create a root system, even when there is only one orphan */
	map<string,MecSystem*,less<string> > mChilds;
	/** Collect orphan systems */
	for(int j=0; j<vsaux.size(); j++)
	{
		mChilds[vsaux[j]->name] = vsaux[j];
	}
	/** Create a new system named "root" with the orphans */
	/** Insert the root system just created */
	mMecSystems["root"] = new MecSystem("root",mChilds, this, dimension);

	rootSys = (mMecSystems["root"]);
	
	if(vsaux.size()==0)     /** No orphan systems, not possible */
	{
		cout << "*ERROR: it is not possible to build a valid system tree in process " << num << "\n";
	}
	/** End of root system set */
}

// ***

void Process::createPoints()
{
	/** Creation of points */
	if (num==rankk) cout << "	Creating points in each system of process " << num << "...\n";
	if(rootSys!=0)	rootSys->createPoints();	/** Recursive from root system */
}

// ***

void Process::setPPointersCompBodies()
{
	/** Pointers to points at components in bodies */
	if(rootSys!=0)	rootSys->setPPointersCompBodies();	/** Recursive from root system */
}

// ***

void Process::setPPointersLoadBodies()
{
	/** Pointers to points at load components in bodies */
	if(rootSys!=0)	rootSys->setPPointersLoadBodies();	/** Recursive from root system */
}

// ***
/*
void Process::setPPointersCompJoints()
{
	// Pointers to points at components in joints
	if(rootSys!=0)	rootSys->setPPointersCompJoints();	// Recursive from root system
}

// ***

void Process::setPPointersLoadJoints()
{
	// Pointers to points at load components in joints
	if(rootSys!=0)	rootSys->setPPointersLoadJoints();	// Recursive from root system
}
*/
// ***

void Process::setPPointersConsJoints()
{
	/** Pointers to points at load components in joints */
	if(rootSys!=0)	rootSys->setPPointersConsJoints();	/** Recursive from root system */
}

// ***

void Process::setPPointSurface()
{
	/** Pointers to points at load components in joints */
	if(rootSys!=0)	rootSys->setPPointSurface();	/** Recursive from root system */
}

// ***

void Process::setRelProc()
{
	map<string, Point, less<string> >::iterator i;
//JC	map<string, Point*, less<string> >::iterator i;
	string nameSys;
	int procNumber;
	
	relatedProc.set(num,1);			/** First, put 1 at the local process */
	/** Loop on the points inside the root system */
	for(i=(rootSys->mPoints).begin(); i!=(rootSys->mPoints).end(); i++)
	{
		nameSys = getStringHead((i->second).name);
//JC		nameSys = getStringHead(i->second->name);
		procNumber = mSysProc[nameSys];
		relatedProc.set(procNumber,1);
		(getPTProc(nameSys)->relatedProc).set(num,1);		/** Put 1 in the other system */
	}
//	cout << "RelatedProc = " << relatedProc << endl;
}

// ***

vector<Point*> Process::getPTPoints(vector<string>& pointNames)
{
	// Returns a vector of local pointers from a vector of point names
	int i;
//	Point* aux;
	vector<Point*> vP;
//	string namesys;
	
	for(i=0; i<pointNames.size(); i++)
	{
		vP.push_back(getPTPoint(pointNames[i]));
		/*
		namesys=getStringHead(pointNames[i]);

		// Check if the point is local
		if(checkIfSystem(namesys)==true)
		{
			// Check if the name changed
			if ( mMecSystems[namesys]->checkIfPointChanged(pointNames[i])==true )
			{
				// Get the new name
//				if(rankk==num) cout << "En el prpceso " << num << " El punto " << pointNames[i] << " ha cambiado de nombre " << endl;
				pointNames[i] = mMecSystems[namesys]->mChangedPoints[pointNames[i]];
//				if(rankk==num) cout << "El nuevo nombre es " << pointNames[i] << endl;
			}
		}

		// Check if the point is in the process
		if(checkIfSystem(namesys)==true)
		{
			aux = &(mMecSystems[namesys]->mPoints[pointNames[i]]);
		}
		// Look at the root system
		else if(checkIfSystem(namesys)==false)
		{
			aux = &((rootSys)->mPoints[pointNames[i]]);
		}
		vP.push_back(aux);
		 */
	}
	return(vP);
}

// ***

Point* Process::getPTPoint(const string& pointName)
{

	// Returns a local pointer from a single point name
	Point* aux;
	string namesys;
	string auxs;
	
//	if(rankk==num) cout << "Buscando puntero de punto " << pointName << endl;
	namesys=getStringHead(pointName);
	//		cout << "El punto es " << pointNames[i] << endl;
	//		cout << "Namesys es " << namesys << endl;
	
	auxs = pointName;
	
	// Check if the point is local
	if(checkIfSystem(namesys)==true)
	{
		// Check if the name changed
		if ( mMecSystems[namesys]->checkIfPointChanged(pointName)==true )
		{
			// Get the new name
//			if(rankk==num) cout << "En el prpceso " << num << " El punto " << pointName << " ha cambiado de nombre " << endl;
			auxs = mMecSystems[namesys]->mChangedPoints[pointName];
//			if(rankk==num) cout << "El nuevo nombre es " << auxs << endl;
		}
	}
	
	namesys=getStringHead(auxs);
	
	// Check if the point is in the process
	if(checkIfSystem(namesys)==true)
	{
		//						cout << "Es local " <<  mMecSystems[namesys]->mPoints[pointNames[i]].name << endl;
		aux = &(mMecSystems[namesys]->mPoints[auxs]);
	}
	// Look at the root system
	else if(checkIfSystem(namesys)==false)
	{
//		if(rankk==num) cout << "He encontrado en el proceso " << num << " un punto remoto, el " << auxs << endl;
		aux = &((rootSys)->mPoints[auxs]);
	}
	
// -----------
/* 
	namesys=getStringHead(pointName);

	if(checkIfSystem(namesys)==false)
	{
		aux = &((rootSys)->mPoints[pointName]);
	}
	else
	{
		aux = &(mMecSystems[namesys]->mPoints[pointName]);
	}
*/
	return(aux);
}

// ***

void Process::createInterfaces(const int& integratorType)
{

	map<string, int, less<string> >::iterator k;
	map<string, Point, less<string> >::iterator r;
	int i;
	if( relatedProc.count()>1 )
	{
		for(i=0; i<20; i++)
		{
			if( relatedProc.test(i)==true && i!=num )
			{
				vInterfaces.push_back( *(new Interface(num,i,integratorType,this)) );
				//JCarma vInterfaces.push_back( *(new InterfaceArma(num,i,integratorType,this)) );
				//JCeigen vInterfaces.push_back( *(new InterfaceEigen(num,i,integratorType,this)) );
				if (num==rankk) cout << "		Interface between process " << num << " and process " << i << " created \n";
			}
		}
	}
	numInterfaces = vInterfaces.size();
}

// ***

void Process::relateInterfaces()
{
	if( vInterfaces.size()>0)
	{
		int i;
		int	j;
		int k;
		vector<Interface*> vI;
		//JCarma vector<InterfaceArma*> vI;
		//JCeigen vector<InterfaceEigen*> vI;
		vI.resize(numInterfaces-1);

		for(i=0; i<numInterfaces; i++)	(vInterfaces[i]).setRemoteInterface();
		for(i=0; i<numInterfaces; i++)	(vInterfaces[i]).synchronize();
		for(i=0; i<numInterfaces; i++)	(vInterfaces[i]).getPointers();
		
		// Send, to each interface, a vector of pointers to the other interfaces of the local process
		for(i=0; i<numInterfaces; i++)	
		{
			k=0;
			for(j=0; j<numInterfaces; j++)
			{
				if(j!=i)
				{
					vI[k] = &vInterfaces[j];
					k++;
				}
			}
			vInterfaces[i].setOtherInter(vI);
		}
	}
}

// ***

void Process::calcSizeInterfaces(const int& integratorType, bool& implicitFlag)
{
	int i;

	if (implicitFlag==true)	// Implicit problem
	{
		int j;
		int cont;
		int contLG = -1;		

		if (integratorType==2) symmFlag = false;	// Non-symmetric matrices for conserving integration
	
		globalLocalI.resize(numDofs);
		globalLocalB.resize(numDofs);
	
		for(i=0; i<numDofs; ++i)
		{
			globalLocalI[i] = -1;
			globalLocalB[i] = -1;
		}
		
		for(i=0; i<numInterfaces; i++)	(vInterfaces[i]).calcSize();

		numDofsI = numDofs;
	
		if(numInterfaces!=0)
		{
			for( i=0; i<numDofs; ++i) 
			{
				cont = 0;
				for(j=0; j<numInterfaces; j++)
				{
					if ( vInterfaces[j].checkIfDof(i)==true )
					{
						cont++;
					}
				}
				if(cont==0)
				{
					localGlobI.push_back(i);
					contLG+=1;
					globalLocalI[i] = contLG;
				}
				else	localGlobB.push_back(i);
			}
			numDofsI = localGlobI.size();
			numDofsB = localGlobB.size();
		
			for (i=0; i<numDofsB; ++i)
			{
				if( globalLocalI[ localGlobB[i] ] ==-1)	globalLocalB[ localGlobB[i] ] = i;
			}
		
			//JC
			Kbb.resize(numDofsB);
			for (i=0; i<numDofsB; ++i)	Kbb[i].resize(numDofsB);
		
			// B is (-Kib | Ri)
			// B is transposed
			B.resize(numDofsB+1);
			for (i=0; i<numDofsB+1; ++i)	B[i].resize(numDofsI);
			//
			// Kib
			Kib.resize(numDofsI);
			for (i=0; i<numDofsI; ++i)	Kib[i].resize(numDofsB);
			//
			Ri.resize(numDofsI);
			Rb.resize(numDofsB);
			if (symmFlag==false)
			{
				Kbi.resize(numDofsB);
				for (i=0; i<numDofsB; ++i)	Kbi[i].resize(numDofsI);
			}
			//JC
		}
		/*
		 if (num==1)
		 {
		 cout << "globalLocalI is " << endl;
		 printSTLVector(globalLocalI);

		 cout << "globalLocalB is " << endl;
		 printSTLVector(globalLocalB);
		 }
		 */
	}
	else	// Explicit
	{
		for(i=0; i<numInterfaces; i++)	(vInterfaces[i]).calcSize();
	}
}

// ***

void Process::setSizeInterfaces(bool& implicitFlag)
{
	int i;

	if (implicitFlag==true)
	{
//		#pragma omp parallel for private (i)
		for(i=0; i<numInterfaces; i++)	(vInterfaces[i]).setSize();
	}
	else 
	{
//		#pragma omp parallel for private (i)
		for(i=0; i<numInterfaces; i++)	(vInterfaces[i]).setSizeExp();
	}
}

// ***

void Process::addRemoteSharedPoint(string pointName, int numProc)
{
	string sysName = getStringHead(pointName);
	int hostProc = mSysProc[sysName];
	
	(vPTProc[hostProc])->addLocalSharedPoint(pointName,numProc);
}

// ***

void Process::addLocalSharedPoint(string pointName, int numProc)
{
	localSharedPoints.insert( pair<string,int>(pointName,numProc) );
}

// ***

void Process::setDofsBodies()
{	
		rootSys->setDofsBodies();			/** Recursive from root system */
}

// ***

void Process::setDofsRoots()
{	
		// Check for points without dofs, located at root system */
		// These are points that are originally from a remote process, and do not belong to any local body
		// In this case, dofs have not been counted. This is tipically the case of joint constraints, that do not add dofs
		// The number of dofs should be copied from the point at the remote process
		rootSys->setDofsRoot();
}

// ***

void Process::setDofsJoints()
{	
		rootSys->setDofsJoints();			/** Recursive from root system */
}

// ***

void Process::setGlobalDofNumbers()
{
	// Recursive from root system
	numDofs=0;

	// Comment one of the options
/*	
	// Option 1: Rotational dofs at the end of vector of unknowns q
	int numTDofs=0;
	int numRDofs=0;
	if(rootSys!=0)	numTDofs = rootSys->setTransNumbers(0);		// Set global numbers for translational dofs
	if(rootSys!=0)	numRDofs = rootSys->setRotNumbers(0);		// Set global numbers for rotational dofs
	numDofs = numTDofs + numRDofs;
	cout << "	-->	Process " << num << " has " << numTDofs << " Tdofs and " << numRDofs << " Rdofs \n";
*/
	// Option 2: Translational and Rotational dofs mixed

	if(rootSys!=0)	numDofs = rootSys->setTRNumbers(0);			// Set global numbers for translational and rotational dofs
	if (num==rankk) cout << "	--> Process " << num << " has " << numDofs << " dofs \n";
	
	// Creates an vector of position and velocities unknowns q, dq with the proper size
	q.resize(numDofs);
	dq.resize(numDofs);
	qIni.resize(numDofs);
	
	//if (dynFlag!=0) setInertiaDofNumbers();
	setInertiaDofNumbers();
}

// ***

void Process::setLocalGlob()
{
	// Recursive from root system
	if(rootSys!=0)	rootSys->setLocalGlob();
//	if (projFlag==true)	setConstraintNumbers();
}

// ***

void Process::sendLoadToBody(string lType, string sName, string bName, ifstream &inn, const vector<double> & vTime) // edit
{
	mMecSystems[sName]->mBodies[bName]->setLoad(lType, inn, vTime, &mMecSystems[sName]->pointSets);
}

// ***

void Process::initializeGlobalDofsBodies()
{
	double tt=0.0;
	rootSys->initializeGlobalDofsBodies(q,dq,tt);
	
	/*
	if(rootSys!=0)
	{
		if(dynFlag==0) rootSys->initializeGlobalDofsBodies(q);
		else rootSys->initializeGlobalDofsBodies(q,dq,tt);
	}
	 */
}

// ***

void Process::initializeGlobalDofsJoints()
{
	double tt=0.0;
	rootSys->initializeGlobalDofsJoints(q,dq,tt);
	
	/*
	if(rootSys!=0)
	{
		if(dynFlag==0) rootSys->initializeGlobalDofsJoints(q);
		else rootSys->initializeGlobalDofsJoints(q,dq,tt);
	}
	*/
}	

// ***

void Process::synchronizeSharedDofs()
{
	
	multimap<string, int, less<string> > ::iterator i;
	vector<int> dofs;
	vector<int> remoteDofs;
	vector<double> valueDofs;
	int j;
	Process* remoteProcess;
	
	// Going through the localSharedPoints
	for(i=localSharedPoints.begin(); i!=localSharedPoints.end(); ++i)
	{
		dofs = getPTPoint( i->first )->getGlobalNumbers();	// Get global local dof numbers
		for (j=0; j<dofs.size(); ++j)
		{
			valueDofs.push_back( q[ dofs[j] ] );				// Get local dofs values
		}
		remoteProcess = getPTProc(i->second);
		remoteDofs = remoteProcess->getPTPoint(i->first)->getGlobalNumbers();	// Get remote dof numbers
		remoteProcess->setDofValues(remoteDofs, valueDofs);						// Set remote dofs values
		valueDofs.clear();
	}
}


// ***

void Process::setVectorQ(vector<double> & q_in)
{
	q = q_in;
}

// ***

void Process::setMaterials(const map<string, Material*, less<string> >& mMat)
{
	mMaterials = mMat;
}

// ***

void Process::setMatAtComp()
{
	// Recursive from root system
	if(rootSys!=0)	rootSys->setMatAtComp();
}

// ***

void Process::formIsolatedResidue(const double& time, const int& modeFlag, const int& modeCalc, const int& integratorType, const double& DeltaT)
{
	setZero(res);
	eneResMax = 0.0;
	int i;
	
//	cout << "ddq en formIsolatedResidue ANTES es " << endl;
//	printSTLVector(ddq);
	
//	cout << "dq en formIsolatedResidue ANTES es " << endl;
//	printSTLVector(dq);


	if(dynFlag==0)	// Statics. Always relative residual energy
	{
		rootSys->fillResidual(res, q, dq, time, modeFlag, grav, gravDir, DeltaT);
	} else	// Dynamics
	{
		if(modeCalc==0)	// Relative residual energy
		{
//JC			rootSys->fillResidual(res, q, dq, ddq, time, modeFlag, grav, gravDir, integratorType, qConv, dqConv, DeltaT);
				rootSys->fillResidual(res, q, dq, ddq, time, modeFlag, grav, gravDir, integratorType, DeltaT);
		}
		else if(modeCalc==1)	// Balanced relative residual energy
								// WARNING: need more testing.  Probably does not work well at this time.
		{
			setZero(resFq);
			setZero(resFt);
		
			rootSys->fillResFi(res, q, ddq, modeFlag);
			
//*			eneResFi = fabs(res*Deltaq);
			eneResFi = 0.0;
			for(i=0; i<Deltaq.size(); ++i)
			{
				eneResFi+= res[i]*Deltaq[i];
			}
			eneResFi = fabs(eneResFi);
			
			if(eneResFi>=eneResMax) eneResMax = eneResFi;
		
//			rootSys->fillResFq(resFq, q, dq, modeFlag, integratorType, qConv, dqConv, DeltaT);
			rootSys->fillResFq(resFq, q, dq, modeFlag, integratorType, DeltaT);
			
//*			eneResFq = fabs(resFq*Deltaq);
			eneResFq = 0.0;
			for(i=0; i<Deltaq.size(); ++i)
			{
				eneResFq+= resFq[i]*Deltaq[i];
			}
			eneResFq = fabs(eneResFq);
			
			if(eneResFq>=eneResMax) eneResMax = eneResFq;
		
			rootSys->fillResFt(resFt, time, modeFlag, grav, gravDir);
			
//*			eneResFt = fabs(resFt*Deltaq);
			eneResFt = 0.0;
			for(i=0; i<Deltaq.size(); ++i)
			{
				eneResFt+= resFt[i]*Deltaq[i];
			}
			eneResFt = fabs(eneResFt);
			
			if(eneResFt>=eneResMax) eneResMax = eneResFt;
		
//*			res+= resFq + resFt;
			for(i=0; i<res.size(); ++i)
			{
				res[i]+= resFq[i] + resFt[i];
			}
			
		}
	}

//	if(num==0)
//	{
//		cout << "q en formIsolatedResidue para num " << num << " es " << endl;
//		printSTLVector(q);
//		//cout << "ddq en formIsolatedResidue es " << endl;
//		//printSTLVector(ddq);
//		cout << "res en formIsolatedResidue para num " << num << " es " << endl;
//		printSTLVector(res);
//	}

}

// ***

void Process::formIsolatedTangent(double& time, const int& modeFlag)
{	
	
	Kint->reset();
	setZero(Kbb);
	setZero(Kib);
	setZero(Kbi);

	rootSys->fillTangent(*Kint, Kbb, Kib, Kbi, globalLocalI, globalLocalB, q, dq, time, modeFlag);
	
	int i;
	int j;
	int a;
	int b;
	/*
	if(num==0)
	{
		cout << "Kint en formIsolatedTangent para num " << num << " es " << endl;
		for(i=0; i<numDofsI; ++i)
		{
			for(j=0; j<numDofsI; ++j)				cout << " " << Kint->readElement(i,j);
			cout << endl;
		}
	}	
	*/
	/*
	if (num==0)
	{
		cout << "Kbb is " << endl;
		printSTLMatrix(Kbb);
		cout << "Kib is " << endl;
		printSTLMatrix(Kib);
		cout << "Kbi is " << endl;
		printSTLMatrix(Kbi);
	}
	*/
	
	// Dofs to be removed
	for (i=0; i<vNoDofs.size(); ++i)
	{
		if (vInterfaces.size()!=0)	// More than one processes
		{
			a = globalLocalI[vNoDofs[i]];
			if (a!=-1)
			{
				if ( (time <= vNoDofsEndTimes[i]) || (vNoDofsEndTimes[i] == 0.0) )
				{
					// Loop in rows and columns
					for (j=0; j<numDofs; ++j)
					{
						b = globalLocalI[j];
						if (b!=-1)
						{
							Kint->writeElement(0.0,a,b);
							Kint->writeElement(0.0,b,a);						
						}
					}
					// Loop on diagonal
					Kint->writeElement(1.0,a,a);
				}
			}
		}
		else	// Just one process
		{
			a = vNoDofs[i];
			if (a!=-1)
			{
				if ( (time <= vNoDofsEndTimes[i]) || (vNoDofsEndTimes[i] == 0.0) )
				{
					// Loop in rows and columns
					for (j=0; j<numDofs; ++j)
					{
						Kint->writeElement(0.0,a,j);
						Kint->writeElement(0.0,j,a);						
					}
					// Loop on diagonal
					Kint->writeElement(1.0,a,a);
				}
			}
		}
	}
	
	// Imposed dofs
	for (i=0; i<vImposedDofs.size(); ++i)
	{
		if (numInterfaces!=0)	// More than one processes
		{
			a = globalLocalI[vImposedDofs[i]];
			if (a!=-1)
			{
				// Check if the dof is imposed at this time step 
				if ( mHistories[vImposedDofsHis[i]]->value(time) != -9999.9999)
				{					
					// Loop in rows and columns
					for (j=0; j<numDofs; ++j)
					{
						b = globalLocalI[j];
						if (b!=-1)
						{
							Kint->writeElement(0.0,a,b);
						}
					}
					// Loop on diagonal
					Kint->writeElement(1.0,a,a);
				}
			}
		}
		else	// Just one process
		{
			a = vImposedDofs[i];
			if (a!=-1)
			{
				// Check if the dof is imposed at this time step 
				if ( mHistories[vImposedDofsHis[i]]->value(time) != -9999.9999)
				{					
					// Loop in rows and columns
					for (j=0; j<numDofs; ++j)
					{
						Kint->writeElement(0.0,a,j);
					}
					// Loop on diagonal
					Kint->writeElement(1.0,a,a);
				}			
			}
		}
	}		
}

// ***

void Process::formIsolatedTangent(double& time, double &partial_dq, double &partial_ddq, const int& modeFlag, const int& integratorType, const double& DeltaT)
{
	// Dynamics
	
	Kint->reset();
	setZero(Kbb);
	setZero(Kib);
	setZero(Kbi);

	rootSys->fillTangent(*Kint, Kbb, Kib, Kbi, globalLocalI, globalLocalB, q, dq, ddq, partial_dq, partial_ddq, time, modeFlag, integratorType, DeltaT);
	
	int i;
	int j;
	int a;
	int b;
	/*
	if(num==0)
	{
		cout << "Kint en formIsolatedTangent para num " << num << " es " << endl;
		for(i=0; i<numDofsI; ++i)
		{
			for(j=0; j<numDofsI; ++j)				cout << " " << Kint->readElement(i,j);
			cout << endl;
		}
	}	
	*/
	// Dofs to be removed
	for (i=0; i<vNoDofs.size(); ++i)
	{
		if (numInterfaces!=0)	// More than one processes
		{
			a = globalLocalI[vNoDofs[i]];
			if (a!=-1)
			{
				if ( (time <= vNoDofsEndTimes[i]) || (vNoDofsEndTimes[i] == 0.0) )
				{
					// Loop in rows and columns
					for (j=0; j<numDofs; ++j)
					{
						b = globalLocalI[j];
						if (b!=-1)
						{
							Kint->writeElement(0.0,a,b);
							Kint->writeElement(0.0,b,a);						
						}
					}
					// Loop on diagonal
					Kint->writeElement(1.0,a,a);
				}
			}
		}
		else	// Just one process
		{
			a = vNoDofs[i];
			if (a!=-1)
			{
				if ( (time <= vNoDofsEndTimes[i]) || (vNoDofsEndTimes[i] == 0.0) )
				{
					// Loop in rows and columns
					for (j=0; j<numDofs; ++j)
					{
						Kint->writeElement(0.0,a,j);
						Kint->writeElement(0.0,j,a);						
					}
					// Loop on diagonal
					Kint->writeElement(1.0,a,a);
				}
			}
		}
	}
	// Imposed dofs
	for (i=0; i<vImposedDofs.size(); ++i)
	{
		if (numInterfaces!=0)	// More than one processes
		{
			a = globalLocalI[vImposedDofs[i]];
			if (a!=-1)
			{
				// Check if the dof is imposed at this time step 
				if ( mHistories[vImposedDofsHis[i]]->value(time) != -9999.9999)
				{
					// Loop in rows and columns
					for (j=0; j<numDofs; ++j)
					{
						b = globalLocalI[j];
						if (b!=-1)
						{
							Kint->writeElement(0.0,a,b);
						}
					}
					// Loop on diagonal
					Kint->writeElement(1.0,a,a);
				
				}
			}
		}
		else	// Just one process
		{
			a = vImposedDofs[i];
			if (a!=-1)
			{
				// Check if the dof is imposed at this time step 
				if ( mHistories[vImposedDofsHis[i]]->value(time) != -9999.9999)
				{
					// Loop in rows and columns
					for (j=0; j<numDofs; ++j)
					{
						Kint->writeElement(0.0,a,j);
					}
					// Loop on diagonal
					Kint->writeElement(1.0,a,a);
				}
			}
		}
	}			
}


// ***

void Process::resetInterInt()
{
	int i;
	// Loop on interfaces, reset internal contributions
//	#pragma omp parallel for private (i)
	for(i=0; i<numInterfaces; i++)	vInterfaces[i].resetInterInt();
}

// ***

void Process::factorizeInternalMatrix(int facMode)
{
	Kint->factorize(facMode);
}

// ***

void Process::setInterExt(bool b)
{
	int i;

//#pragma omp parallel for private (i)
	for(i=0; i<numInterfaces; i++)
	{
		if (b==true)
		{
			vInterfaces[i].initializeKR(Kbb, Kib, Kbi, globalLocalI, globalLocalB, res, *Kint, matType);
		}
		vInterfaces[i].setInterExt(Kbb, Kib, Kbi, globalLocalI, globalLocalB, res, *Kint, matType, b);
	}
//	cout << "External contributions of interfaces of process " << num << " computed" << endl;
}

// ***

void Process::setInterInt()
{
	int i;
	
	// WARNING setInterInt can not be parallelilzed with openmp
	// Loop on interfaces, compute internal contributions
	for(i=0; i<numInterfaces; i++)
	{
		vInterfaces[i].setInterInt();
	}
//	cout << "Internal contributions at interfaces of process " << num << " already set" << endl;
}

// ***

void Process::formResidue(const int& modeCalc)
{
	int i;
	if(dynFlag==0)	// Statics. Always relative residual energy
	{
//		#pragma omp parallel for private (i)
		for(i=0; i<numInterfaces; i++)	vInterfaces[i].addRes(res);
	} else	// Dynamics
	{
		if(modeCalc==0)	// Relative residual energy
		{
//			#pragma omp parallel for private (i)
			for(i=0; i<numInterfaces; i++)	vInterfaces[i].addRes(res);
		}
		else if(modeCalc==1)	// Balanced relative residual energy
		{			
			if(numInterfaces!=0)
			{
				setZero(resInter);
				for(i=0; i<numInterfaces; i++)	vInterfaces[i].addRes(resInter);
				
//*				eneResInter = fabs(resInter*Deltaq);
				eneResInter = 0.0;
				for(i=0; i<Deltaq.size(); ++i)
				{
					eneResInter+= resInter[i]*Deltaq[i];
				}
				eneResInter = fabs(eneResInter);
				
				if(eneResInter>=eneResMax) eneResMax = eneResInter;
				
//*				res+= resInter;
				for(i=0; i<res.size(); ++i)
				{
					res[i]+= resInter[i];
				}				
			}
//*			eneRes = fabs(res*Deltaq);
			eneRes = 0.0;
			for(i=0; i<Deltaq.size(); ++i)
			{
				eneRes+= res[i]*Deltaq[i];
			}
			eneRes = fabs(eneRes);
		}
	}
/*	
	if(num==0)
	{
		cout << "res en formResidue es " << res.norm2() << endl;
	}
*/
	/*
	cout << "Residual is " << endl;
	printSTLVector(res);
	*/
	// Fill Rb
	//JC
	if(numInterfaces!=0)
	{
		// Ri is the last column of matrix B, filled in formTangent, but still I need it
		for(i=0; i<numDofsI; ++i)	Ri[i] = res[localGlobI[i]];
		for(i=0; i<numDofsB; ++i)	Rb[i] = res[localGlobB[i]];
	}
}

// ***

void Process::formTangent()
{
	int i;
	int j;
	
	// For each interface, add the internal residues to res
//	#pragma omp parallel for private (i)
	for(i=0; i<numInterfaces; i++)
	{
//		vInterfaces[i].addTang(*tang,matType);
		vInterfaces[i].addTang(Kbb, Kib, Kbi, globalLocalI, globalLocalB, matType);
	}

	if(numInterfaces!=0)
	{
		// B = (-Kib | Ri)
		// B is transposed
		for(i=0; i<numDofsI; ++i)
		{
			for (j=0; j<numDofsB; ++j)	B[j][i] = - Kib[i][j];
			B[numDofsB][i] = res[localGlobI[i]];
		}
	}
	/*
	if (num==0)
	{
		cout << "Kbb = " << endl;
		printSTLMatrix(Kbb);
		cout << "Kib = " << endl;
		printSTLMatrix(Kib);
		cout << "Kbi = " << endl;
		printSTLMatrix(Kbi);
	}
	*/
	/*
	cout << "Tangent is " << endl;
	for (i=0; i<numDofs; ++i)
	{
		for (j=0; j<numDofs; ++j)
		{
			cout << Kint->readElement(i,j) << "  ";
		}
		cout << endl;
	}
    */
}

// ***

void Process::myEvaluationImplicit	(const double& time, const int& modeCalc, const bool& calcM)
{
	
	// Only for implicit integrators
	

	if(modeCalc==1 )
	{
		if (Deltaq.size() == 0){Deltaq.resize(numDofs);}
	}
	
	// ddq has initially the internal+external forces
	if (ddq.size() == 0)	ddq.resize(numDofs);

	sparseMatrix* Mfac;
	if (calcM==true) 
	{
		if (matType==1)
		{
			Mfac = new sparseMatrixSLU();
		}
		if (matType==2)
		{
			Mfac = new sparseMatrixUMF();
		}
		setPatternInertia(Mfac);
		rootSys->fillMassMatrix(*Mfac,q,time);
		/*
		cout << "NumDofs is " << numDofs <<  " " << Mfac->n << endl;
		cout << "Mass matrix is " << endl;
		for (int i=0; i<Mfac->n; ++i)
		{
			for (int j=0; j<Mfac->n; ++j)
			{
				cout << Mfac->readElement(i, j) << " ";
			}
			cout << endl;
		}
		*/
		if (projFlag==true)	Mfac->copyMatrix(M);
		if (autoProjFlag!=0)	Mfac->copyMatrix(MfacProj);

		// Remember that M is numDofsInertia x numDofsInertia
		Mfac->factorize(0);	// 0 means factorize from scratch

		// This must be optimized, such that only one factorization is performed. **** TO BE DONE
		if (autoProjFlag!=0)	MfacProj->factorize(0);	
	}

	setZero(ddq);
	rootSys->fillForce(ddq,q,dq,time,grav,gravDir);

	if(numDofsInertia!=numDofs)
	{
		// cout << "numDofs = " << numDofs << " ; numDofsInertia = " << numDofsInertia << endl;
		// Creating a new vector ddq that matches the dimensions of M
		// Numbers in numNoInertia are ordered
		vector<double> ddqAux;
//		cout << "ddq = " << endl;
//		printSTLVector(ddq);setSpars
//		cout << "ddq.size = " << ddq.size() << endl;
		removeFromVector(ddq, numNoInertia, ddqAux);
//				cout << "ddqAux.size = " << ddqAux.size() << endl;
//		cout << "ddqAux antes = " << endl;
//		printSTLVector(ddqAux);
		Mfac->subsSolve(ddqAux);
//		M->subsSolve(ddqAux);
//		cout << "ddqAux despues = " << endl;
//		printSTLVector(ddqAux);
		restoreFromVector(ddq, numNoInertia, ddqAux);
//		cout << "ddq despues = " << endl;
//		printSTLVector(ddq);
	}
	else
	{
		Mfac->subsSolve(ddq);
	}
	
	/*
	int i;
	int h=0.0;
	// Initial acceleration of imposed dofs
	for (i=0; i<vImposedDofs.size(); ++i)
	{
		if (h!=0.0)
		{
			ddq[vImposedDofs[i]] =  1.0e7;
			//cout << "dq[] = " << dq[vImposedDofs[i]] << endl;
		}
	}	
	*/
}

// ***

void Process::myEvaluationExplicit	(const double& time, const bool& calcM, const bool& lumped)
{
	// Only for explicit integrators.
	
	// calcM = true means that a consistent M has to be computed
	int i;
	
	if (ddq.size() == 0)	ddq.resize(numDofs);

	if (calcM==true) 
	{
		if (lumped==false)	// Consistent mass matrix
		{
			if (matType==1)
			{
				MfacExp = new sparseMatrixSLU();
			}
			if (matType==2)
			{
				MfacExp = new sparseMatrixUMF();
			}
			setPatternInertia(MfacExp);
			rootSys->fillMassMatrix(*MfacExp,q,time);
		
			if (projFlag==true)	MfacExp->copyMatrix(M);
			if (autoProjFlag!=0)	MfacExp->copyMatrix(MfacProj);
		
			// Remember that M is numDofsInertia x numDofsInertia
			MfacExp->factorize(0);	// 0 means factorize from scratch
		
			// This must be optimized, such that only one factorization is performed. **** TO BE DONE
			if (autoProjFlag!=0)	MfacProj->factorize(0);
		}
		//else
		//{
			// Already created
			//vLumpedMass.resize(numDofsInertia);
			//rootSys->fillMassMatrix(vLumpedMass,q,time);
		//}
	}
	
	/*
	cout << "Mass matrix lumped en proceso " << num << endl;
	printSTLVector(vLumpedMass);
	*/
	
	setZero(ddq);
	
	for (i=0; i<vNoDofs.size(); ++i)
	{
		if ( (time <= vNoDofsEndTimes[i]) || (vNoDofsEndTimes[i] == 0.0) )
		{
			q[vNoDofs[i]] = qConv[vNoDofs[i]];
			dq[vNoDofs[i]] = 0.0;
		}
	}
	
	rootSys->fillForce(ddq,q,dq,time,grav,gravDir);
	/*
	cout << "***  fuerzas aisladas en proceso " << num << endl;
	printSTLVector(ddq);
	*/
	
//	#pragma omp parallel for private (i)
	for(i=0; i<numInterfaces; i++)	vInterfaces[i].initializeRExp(ddq);

	// WARNING setInterIntExp can not be parallelilzed with openmp
	for(i=0; i<numInterfaces; i++)	vInterfaces[i].setInterIntExp();
	
//	#pragma omp parallel for private (i)
	for(i=0; i<numInterfaces; i++)	vInterfaces[i].addResExp(ddq);
	
	/*
	for(i=0; i<numInterfaces; i++)
	{
		vInterfaces[i].initializeRExp(ddq);
		vInterfaces[i].setInterIntExp();
		vInterfaces[i].addResExp(ddq);
	}
	*/
	/*
	cout << "***  fuerzas con interfaces en proceso " << num << endl;
	printSTLVector(ddq);
	*/
	
	// Just in case there are dofs without inertia (e.g. temperature)
	if(numDofsInertia!=numDofs)
	{
		// Creating a new vector ddq that matches the dimensions of M
		// Numbers in numNoInertia are ordered
		vector<double> ddqAux;
		removeFromVector(ddq, numNoInertia, ddqAux);
		if (lumped==false)	MfacExp->subsSolve(ddqAux);	// Consistent
		else	//Lumped
		{
			int i;
			for (i=0; i<numDofsInertia; ++i)	ddqAux[i] /= vLumpedMass[i]; 
		}
		restoreFromVector(ddq, numNoInertia, ddqAux);
	}
	else
	{
		if (lumped==false)	MfacExp->subsSolve(ddq);
		else	//Lumped
		{
			for (i=0; i<numDofsInertia; ++i)	ddq[i] /= vLumpedMass[i]; 
		}		
	}
	
	/*
	cout << "***  ddq (aceleracion) en proceso " << num << endl;
	printSTLVector(ddq);
	*/
}

// ***

void Process::solveUpdate(vector<double>& v, int facMode, double& tim)
{
	int i;
	int a;
	double auxd;

	v = res;
	/*
	if (num==0)
	{
		cout << "residue in process " << num << " is " << endl;
		printSTLVector(v);

		cout << "tang in process " << num << " is " << endl;
		printSparseMatrix(*Kint);
	}
    int kki=8;
    cout << (Kint->readElement(kki,kki))*1e4 << endl;
	*/
	
	// Set the value of imposed dofs at the residual
	// Assuming the imposed dofs are internal

	for (i=0; i<vImposedDofs.size(); ++i)
	{
		if (vInterfaces.size()!=0)	// More than one processes
		{
			a = globalLocalI[vImposedDofs[i]];
		}
		else
		{
			a = vImposedDofs[i];
		}
		if (a!=-1)
		{
			auxd = mHistories[vImposedDofsHis[i]]->value(tim);
			if ( auxd!= -9999.9999)
			{
				v[a] = -(qIni[vImposedDofs[i]] + auxd - q[vImposedDofs[i]]);
			}
		}
		/*	
		else	// Just one process
		{
			a = vImposedDofs[i];
			if (a!=-1)
			{
				//cout << vImposedDofsHis[i] << endl;
				//cout <<  mHistories[ vImposedDofsHis[i] ]->value(tim) << endl;
				//v[a] = -mHistories[ vImposedDofsHis[i] ]->value(tim);
				v[a] = -(qIni[vImposedDofs[i]] + mHistories[ vImposedDofsHis[i] ]->value(tim) - q[vImposedDofs[i]]);
				
				//cout << "Para time = " << tim << " y dof " << a << ", qIni es " <<  qIni[vImposedDofs[i]] << endl;
				//cout << "Para time = " << tim << " y dof " << a << ", se lee en el fichero " << mHistories[ vImposedDofsHis[i] ]->value(tim) << endl;
				//cout << "Para time = " << tim << " y dof " << a << ", el dof q vale " <<  q[vImposedDofs[i]] << endl;
				//cout << "Para time = " << tim << " y dof " << a << ", el incremento debe ser " <<  v[a] << endl;
			}
		}
		*/
	}	
	
	
	if (numInterfaces!=0)
	{
		Kint->subsSolveTrans(B);	// B is now (Q|X) transposed
		// Kbb is overwritten by k = Kbb + Kbi*Q
		int i, j, k;
		for (i=0; i<numDofsB; ++i)
		{
			for (j=0; j<numDofsB; ++j)
			{
				for (k=0; k<numDofsI; ++k)
				{
					if (symmFlag==false)	Kbb[i][j] += Kbi[i][k]*B[j][k];
					else Kbb[i][j] += Kib[k][i]*B[j][k];
				}
			}
		}
		// Rb is overwritten by r = Rb-Kbi*X
		// X is the last row of B (remember that B is transposed)
		for (i=0; i<numDofsB; ++i)
		{
			for (j=0; j<numDofsI; ++j)
			{
				if (symmFlag==false)	Rb[i]-= Kbi[i][j]*B[numDofsB][j];
				else Rb[i]-= Kib[j][i]*B[numDofsB][j];
			}
		}
		lu(Kbb);
		forSub(Kbb,Rb);
		backSub(Kbb,Rb);	// Rb contains qb
		
		//cout << "qb = " << endl;
		//printSTLVector(Rb);
		
		// Ri is overwritten by Ri-Kib*qb
		for (i=0; i<numDofsI; ++i)
		{
			for (j=0; j<numDofsB; ++j)
			{
				Ri[i]-= Kib[i][j]*Rb[j]; 
			}
		}
		
		Kint->subsSolve(Ri);	// Ri contains qi
		
		//cout << "qi = " << endl;
		//printSTLVector(Ri);
		
		// Re-construct v
		for (i=0; i<numDofsI; ++i)	v[localGlobI[i]] = Ri[i];
		for (i=0; i<numDofsB; ++i)	v[localGlobB[i]] = Rb[i];		
	}
	else 
	{
		Kint->subsSolve(v);
	}
	
	/*
	if(num==0)
	{
	cout << "delta q in process " << num << " is " << endl;
	printSTLVector(v);


	cout << "old q in process " << num << " is " << endl;
	printSTLVector(q);
	}
	*/
	
	// Put 0 increment in the deleted dofs, and 0 res
	for (i=0; i<vNoDofs.size(); ++i)
	{
		//cout << "Poniendo 0  en vez de " << v[vNoDofs[i]] << " en el lugar " << vNoDofs[i] << " de v " << endl;
		if ( (tim <= vNoDofsEndTimes[i]) || (vNoDofsEndTimes[i] == 0.0) )
		{
			v[vNoDofs[i]] = 0.0;
			res[vNoDofs[i]] = 0.0;
		}
	}
	
	rootSys->updateDofs(q,v);

	// Set the imposed dofs
	for (i=0; i<vImposedDofs.size(); ++i)
	{
		auxd = mHistories[vImposedDofsHis[i]]->value(tim);
		if ( auxd!= -9999.9999)
		{
			q[vImposedDofs[i]] =  qIni[vImposedDofs[i]] + auxd;
		}
	}
	/*
	cout << "new q in process " << num << " is " << endl;
	printSTLVector(q);
	*/

	if(Deltaq.size()!=0)
	{
		for(i=0; i<q.size(); ++i)	Deltaq[i] = q[i] - qConv[i];
	}
	/*
	arma::mat A = arma::randu<arma::mat>(100,200);
	arma::mat B = arma::randu<arma::mat>(100,200);
	cout << A*arma::trans(B) << endl;
	*/
/*
	mat A = randu<mat>(100,200);
	mat B = randu<mat>(100,200);
	cout << A*trans(B) << endl;
*/	
}

// ***

void Process::reportContent()
{
	cout << "Process " << num << " : \n";
	if(rootSys!=0)	rootSys->reportContent();	/** Recursive from root system */
	else	cout << "*ERROR: Process " << num << " has a no valid system tree \n";
}

// ***

void Process::reportPointsDofs()
{
	rootSys->reportPointsDofs();	/** Recursive from root system */
}

// ***

void Process::reportPointsDofsFile()
{
	rootSys->reportPointsDofsFile();	/** Recursive from root system */
}

// ***

int Process::getNumberInterfaces()
{
	return(numInterfaces); 
}

// ***

void Process::writeToFile(const string& st)
{
	if(num==rankk)
	{
		ou.open(outFileName,ios::app);
		ou << st << endl;
		ou.close();
	}
}

// ***

void Process::writeToFile0(const string& st)
{
	if(num==rankk)
	{
		ou.open(outFileName,ios::app);
		ou << st;
		ou.close();
	}
}

// ***

void Process::writeToFile(const string& st, const char* c)
{
	if(num==rankk)
	{
		ou.open(outFileName,ios::app);
		ou << st << c << endl;
		ou.close();
	}
}

// ***

void Process::writeToFile(const string& st1, const string& st2, const string& st3)
{
	if(num==rankk)
	{
		ou.open(outFileName,ios::app);
		ou << st1 << st2 << st3 << endl;
		ou.close();
	}
}

// ***

void Process::writeToFile(const string& st, const int& i)
{
	if(num==rankk)
	{
		ou.open(outFileName,ios::app);
		ou << st << i << endl;
		ou.close();
	}
}

// ***

void Process::writeToFile(const int& i)
{
	if(num==rankk)
	{
		ou.open(outFileName,ios::app);
		ou << i;
		ou.close();
	}
}

// ***

void Process::writeToFile(const int& i, const string& st)
{
	if(num==rankk)
	{
		ou.open(outFileName,ios::app);
		ou << i << st;
		ou.close();
	}
}

// ***

void Process::writeToFile(const string& st, const int& i, const string& st2, const double& d)
{
	if(num==rankk)
	{
		ou.open(outFileName,ios::app);
		ou << st << i << st2 << d << endl;
		ou.close();
	}
}

// ***

void Process::writeToFile(const vector<double> dd)
{
	if(num==rankk)
	{
		ou.open(outFileName,ios::app);
		ou.close();
	}
}

// ***

void Process::writeToFile(const string& st1, const int& i1, const int& i2, const string& st2)
{
	if(num==rankk)
	{
		ou.open(outFileName,ios::app);
		ou << st1 << i1 << i2 << st2 << endl;
		ou.close();
	}
}

// ***

void Process::writeToFile(const vector< vector<double> >& a)
{
	int i;
	int j;
	int nrows = a.size();
	int ncols = a[0].size();
	
	if(num==rankk)
	{
		ou.open(outFileName,ios::app);
		for (i=0; i<nrows; ++i)
		{
			for (j=0; j<ncols; ++j)
			{
				ou << a[i][j] << "\t";
			}
			ou << "\n";
		}
		ou.close();
	}
}

// ***

void Process::writeToFile(const vector<int>& v)
{
	int i;
	
	if(num==rankk)
	{
		ou.open(outFileName,ios::app);
		for (i=0; i<v.size(); ++i)
		{
			ou << v[i] << "\t";
		}
		ou << "\n";
		ou.close();
	}
}

// ***

void Process::writeToFileTime(const string& st, const double& i)
{
	if(num==rankk)
	{
		ouTime.open(timeInfoFileName,ios::app);
		ouTime << st << i << endl;
		ouTime.close();
	}
}

// ***

void Process::createDofOutFile()
{	
	stringstream stnum;
	string snum; 

	// Create if the output file does not exist
	if ( !strcmp(outDofsFile,"") )
	{
		stnum << num;
		stnum >> snum;
	
		strcpy(outDofsFile,"Dofs-");
		strcat(outDofsFile,title);
		strcat(outDofsFile,"-");
		strcat(outDofsFile,snum.c_str());
	
		stnum.clear();
	
		stnum << rankk;
		stnum >> snum;
		
		strcat(outDofsFile,".dat");
		
		fDofs.open(outDofsFile);
		fDofs.close();
		fDofs.open(outDofsFile,ios::app);
	}
	if ( !strcmp(outVelDofsFile,"") )
	{
		stnum.clear();
		snum.clear();
		
		stnum << rankk;
		stnum >> snum;

		strcpy(outVelDofsFile,"VelDofs");
		strcat(outVelDofsFile,snum.c_str());
		
		stnum.clear();
		
		stnum << rankk;
		stnum >> snum;
		
		strcat(outVelDofsFile,".dat");
		
		fVelDofs.open(outVelDofsFile);		
		fVelDofs.close();
		fVelDofs.open(outVelDofsFile,ios::app);
	}
}

// ***

void Process::createStressOutFile()
{	
	stringstream stnum;
	string snum; 
	
	// Create if the output file does not exist
	if ( !strcmp(outStressFile,"") )
	{
		stnum << num;
		stnum >> snum;
		
		strcpy(outStressFile,"Stress");
		strcat(outStressFile,snum.c_str());
		
		stnum.clear();
		
		stnum << rankk;
		stnum >> snum;
		
		strcat(outStressFile,".dat");
		
		fStress.open(outStressFile);
		fStress.close();
		fStress.open(outStressFile,ios::app);
	}
}

// ***

void Process::fillVDofsOut(const string& pointName)
{
	int i;
	Point* p;
	vector<int> v;
	
//	fDofs.open(outDofsFile,ios::app);
//	fVelDofs.open(outVelDofsFile,ios::app);
		
//	if (rankk==num) cout << "En proceso " << num << " buscando puntero de  " << pointName << endl;
	p = getPTPoint(pointName);
	v = p->getGlobalNumbers();
	
//	cout << "pointName " << pointName << " in process " << num << " with v of size " << v.size() << "\n";
//	printSTLVector(v);
	
	for (i=0; i<v.size(); ++i)
	{
		vDofsOut.push_back(v[i]);
	}

	fDofs << "#" << pointName << endl;
//	fDofs.close();
	
	fVelDofs << "#" << pointName << endl;
//	fVelDofs.close();
}

// ***

void Process::fillVStressOut(const string& pointName)
{
	int i;
	Point* p;
	vector<int> v;
		
	p = getPTPoint(pointName);
	v = p->getGlobalNumbers();
	
	
	for (i=0; i<v.size(); ++i)
	{
		vStressOut.push_back(v[i]);
	}
	
	fStress << "#" << pointName << endl;
}

// ***

void Process::writeFileDofs(const double& time)
{
	int i;
//	fDofs.open(outDofsFile,ios::app);
	
	fDofs << time;
	for (i=0; i<vDofsOut.size(); ++i)
	{
		fDofs << "\t" << q[ vDofsOut[i] ];
	}
	fDofs << "\n";
//	fDofs.close();

	if(dq.size()!=0)
	{
//		fVelDofs.open(outVelDofsFile,ios::app);
		fVelDofs << time;
		for (i=0; i<vDofsOut.size(); ++i)
		{
			fVelDofs << "\t" << dq[ vDofsOut[i] ];
		}
		fVelDofs << "\n";
//		fVelDofs.close();
	}
}

// ***

void Process::writeFileStress(const double& time)
{
	int i;
	int j;
	
	fStress << time;
	for (i=0; i<2; ++i)
	{
		for (j=0; j<vStressOut.size(); ++j)
		{
			fStress << "\t" << FSmooth[i][ vStressOut[j] ];
		}
	}
	for (i=0; i<2; ++i)
	{
		for (j=0; j<vStressOut.size(); ++j)
		{
			fStress << "\t" << FSmoothDerived[i][ vStressOut[j] ];
		}
	}
	
	fStress << "\n";
}

// ***

void Process::writeFileGraph(string &gFormat, string &gPath, const int& tS, bool& plotDispFlag, bool& plotProcFlag, bool& plotStressFlag)
{

	if(!gFormat.empty())
	{
	   fGraph.open(outGraphFile,ios::app);
	   
	   string str;
	   
	   if(!gFormat.compare("vtu"))
	   {
		stringstream stnum;
		string snum;
		stnum << num;
		stnum >> procSS;
		procSS+="-";
	   
		stnum.clear();
		stnum << tS;	// current timestep is added to the name of the graphic outputfile
		stnum >> snum;
	   
		procSS+= snum;	// procSS is a string process-numStep		
	   }
	   if(rootSys!=0)	rootSys->writeFileGraph(fGraph,gFormat,gPath,procSS,plotDispFlag,plotProcFlag, plotStressFlag);	/** Recursive from root system */
	   fGraph.close();
	}
}

// ***

void Process::setDofValues(vector<int>& pos, vector<double>& val)
{
	int i;
	
	// There is no checking that the dimensions match
	for (i=0; i<pos.size(); ++i)
	{
		q[ pos[i] ] = val[i];
	}
}

// ***

void Process::setSparsePattern()
{
		
	cout << "Setting sparse pattern of matrices in process " << num << " ... " << endl;
	
//	cout << "*********** En setSparsePattern \n";
	vector<vector<int> > tempIa;	// Temporal row index
	vector<vector<int> > tempIaI;	// Temporal row index for Kint
	vector<int>::iterator p;
	int i;
	int j;
	int counter;
	int a;
	int b;
//	vector<int> ia;			// Row index (global tangent matrix)
//	vector<int> ja;			// Column index	(global tangent matrix)
//	vector<int> iaInt;		// Row index (internal tangent matrix)
//	vector<int> jaInt;		// Column index (internal tangent matrix)		

	tempIa.resize(numDofs);
	tempIaI.resize(numDofsI);
	
	rootSys->setSparsePattern(tempIa);	/** Recursive from root system */

//	printSTLMatrix(tempIa);

	
	// Add the interface dofs
	// Dofs belonging to the total interface are related	
	for(i=0; i<numDofsB; ++i)
	{
			a = localGlobB[i];
			for(j=0; j<numDofsB; ++j)
			{
				tempIa[ a ].push_back( localGlobB[j] );
			}
	}
	// order and remove duplicates at each row of tempIa

	for (i=0; i<numDofs; ++i)
	{
		sort(tempIa[i].begin(),tempIa[i].end());
		p=unique(tempIa[i].begin(),tempIa[i].end());
		tempIa[i].erase(p,tempIa[i].end());
	}

//	cout << "tempIA global " << endl;
//	printSTLMatrix(tempIa);

	// build column index vector
	// it is assumed that every column has at least one element different from 0
	ja.push_back(1);
	counter = 1;
	for(i=0; i<numDofs; ++i)
	{
		counter += tempIa[i].size();
		ja.push_back( counter );
	}
	
	
	// build row index vector
	for (i=0; i<numDofs; ++i)
	{
		for (j=0; j<tempIa[i].size(); ++j)
		{
			ia.push_back( tempIa[i][j]+1 );
		}
	}
	
	// Define tangent matrix with the computed sparse pattern
	// Sparse pattern of Kint
	if(vInterfaces.size()!=0)
	{
		int A;
		int B;
		for (i=0; i<numDofsI; ++i)
		{
			A = localGlobI[i];
			for (j=0; j<tempIa[A].size(); ++j)
			{
				B = tempIa[A][j];
				b = globalLocalI[B];
				if (b!=-1)	tempIaI[i].push_back(b);
			}
		}
		
		// order and remove duplicates at each row of tempIaI
		for (i=0; i<numDofsI; ++i)
		{
			sort(tempIaI[i].begin(),tempIaI[i].end());
			p=unique(tempIaI[i].begin(),tempIaI[i].end());
			tempIaI[i].erase(p,tempIaI[i].end());
		}
		// build column index vector
		// it is assumed that every column has at least one element different from 0
		jaInt.push_back(1);
		counter = 1;
		for(i=0; i<numDofsI; ++i)
		{
			counter += tempIaI[i].size();
			jaInt.push_back( counter );
		}
		
		
		// build row index vector
		for (i=0; i<numDofsI; ++i)
		{
			for (j=0; j<tempIaI[i].size(); ++j)
			{
				iaInt.push_back( tempIaI[i][j]+1 );
			}
		}
	}
	else 
	{
		iaInt = ia;
		jaInt = ja;
	}

	// Build interface location matrices, to speedup access
	//for(i=0; i<vInterfaces.size(); i++)	(vInterfaces[i]).buildSparseLoc(ia,ja);
	
	
	cout << "Sparse pattern setting done in process " << num << endl;
}

// ***

void Process::setSparsePatternInertia()
{
	//WARNING: It must be called AFTER setSparsePattern
	
	cout << "Setting sparse pattern for inertia of matrices in process " << num << " ... " << endl;

	//cout << "numDofs = " << numDofs << endl;
	//cout << "numDofsInertia = " << numDofsInertia << endl;
	if(numDofs!=numDofsInertia)	 // only if some dofs do not contribute to inertia
	{
		vector<vector<int> > tempIa;	// Temporal row index
		vector<int>::iterator p;
		int i;
		int j;
		int counter;
	
		tempIa.resize(numDofsInertia);			
		rootSys->setSparsePatternInertia(tempIa);	/** Recursive from root system */
			
		// order and remove duplicates at each row of tempIa
		// Each dof without inertia produces an empty row
		for (i=0; i<numDofsInertia; ++i)
		{
			sort(tempIa[i].begin(),tempIa[i].end());
			p=unique(tempIa[i].begin(),tempIa[i].end());
			tempIa[i].erase(p,tempIa[i].end());
		}
			
		// build column index vector
		// it is assumed that every column has at least one element different from 0
		jaInertia.push_back(1);
		counter = 1;
		for(i=0; i<numDofsInertia; ++i)
		{
			if(tempIa[i].size()!=0)
			{
				counter += tempIa[i].size();
				jaInertia.push_back( counter );
			}
		}
			
		// build row index vector
		//			iaInertia.clear();
		for (i=0; i<numDofsInertia; ++i)
		{
			for (j=0; j<tempIa[i].size(); ++j)
			{
				iaInertia.push_back( tempIa[i][j]+1 );
			}
		}
	}
	else
	{
		iaInertia = ia;
		jaInertia = ja;
	}
	/*
	cout << "iaInertia = " << endl;
	printSTLVector(iaInertia);
	cout << "jaInertia = " << endl;
	printSTLVector(jaInertia);
	 */
	cout << "Sparse pattern inertia setting done in process " << num << endl;
}

// ***

void Process::initializeRes(const int& modeCalc, const bool augLagFlag)
{
	res.resize(numDofs);
	if(dynFlag==1 && modeCalc==1)	// Dynamics and balanced residual convergence check
	{
		if(numInterfaces!=0)	resInter.resize(numDofs);
		resFq.resize(numDofs);
		resFt.resize(numDofs);
	}	
	// Initialize residual vectors in components
	rootSys->initializeCompRes(dynFlag,grav,gravDir,augLagFlag);	/** Recursive from root system */
}

// ***

void Process::initializeTang(const int& modeCalc, const bool pf, const int apf, const bool lmf, const bool implicitFlag)
{
	// pf	: Projection flag
	// apf	: Auto-projection flag
	// lmf	: LumpedMassFlag.  true means lumped
	// implicitFlag
	
	int i;
	projFlag = pf;
	autoProjFlag = apf;
	
	//if (matType==1)	tang = new sparseMatrixSLU();
	//if (matType==2)	tang = new sparseMatrixUMF();

	if(dynFlag==0)	// Statics
	{
		if (matType==1)	Kint = new sparseMatrixSLU();
		if (matType==2)	Kint = new sparseMatrixUMF();
		
		setSparsePattern();
		Kint->sparsePattern(iaInt,jaInt);
		Kint->initMatrices();
	}
	else			// Dynamics
	{
		// Internal tangent matrix
		// Only for implicit
		if (implicitFlag==true)
		{
			if (matType==1)	Kint = new sparseMatrixSLU();
			if (matType==2)	Kint = new sparseMatrixUMF();
		}	
		
		if (projFlag==true)	// with velocity projection
		{
			if (lmf==false)	// false means consistent mass matrix
			{
				if (matType==1)
				{
					M = new sparseMatrixSLU();
					pMatrix = new sparseMatrixSLU();
				}
				if (matType==2)
				{
					M = new sparseMatrixUMF();
					pMatrix = new sparseMatrixUMF();		
				}
				if (autoProjFlag!=0)	// Automatical projection (precomputed projection penalty)
				{
					if (matType==1)	MfacProj = new sparseMatrixSLU();
					if (matType==2)	MfacProj = new sparseMatrixUMF();
				}
			}
		}
		// The following parenthesis is not necessary with lumped mass matrix
		if (lmf==false)	// Consistent mass matrix
		{
			if(matType!=0)
			{
				setSparsePattern();
				if (Kint!=0)
				{
					Kint->sparsePattern(iaInt,jaInt);
					Kint->initMatrices();
				}
				setSparsePatternInertia();
				
				if (projFlag==true)
				{
					setPatternInertia(M);
					setPatternInertia(pMatrix);
					projVel.resize(pMatrix->rows());
				}
				if (autoProjFlag!=0)	// Automatical projection (precomputed projection penalty)
				{
					setPatternInertia(MfacProj);
				}
			}
		}		
	}
	
	
			
	// Initialize tangent matrices in components
//	rootSys->initializeCompTang(dynFlag, lmf);	/** Recursive from root system */
	if (implicitFlag==true)	rootSys->initializeCompTang(dynFlag);	/** Recursive from root system */
	/*
	// Prueba de mumps
	vector<int> kki(11);
	vector<int> kkj(6);
	
	kki[0] = 1;
	kki[1] = 4;
	kki[2] = 2;
	kki[3] = 5;
	kki[4] = 1;
	kki[5] = 2;
	kki[6] = 3;
	kki[7] = 4;
	kki[8] = 5;
	kki[9] = 3;
	kki[10] = 5;
	
	kkj[0] = 1;
	kkj[1] = 3;
	kkj[2] = 5;
	kkj[3] = 8;
	kkj[4] = 10;
	kkj[5] = 12;
	
	//printSTLVector(kkj);
	
	sparseMatrix* A;
	A = new sparseMatrixMUMPS();
	A->sparsePattern(kki,kkj);
	A->initMatrices();
	*/
}

// ***

// edit starts here
void Process::createGraphOutFile(const string& gFormat, const string& gPath, const int& tS) 
{

if(!gFormat.empty())
{
	stringstream stnum;
	string snum;
	stnum << num;
	stnum >> procSS;
	
	if(!gFormat.compare("msh"))
	{		
		procSS+="-";
		
		stnum.clear();
		stnum << tS;	// current timestep is added to the name of the graphic outputfile
		stnum >> snum;
		
		procSS+= snum;	// procSS is a string process-numStep
		strcpy(outGraphFile,"Graph");
		strcat(outGraphFile,procSS.c_str());
		
		strcat(outGraphFile,".msh");
	}
	if(!gFormat.compare("vtu") && tS==0)	
	{
        strcpy(outGraphFile, gPath.c_str());
        strcat(outGraphFile,"Graph");
		strcat(outGraphFile,procSS.c_str());
               
		char *dirName = new char [300];
                
        strcpy(dirName, outGraphFile);
               
		procSS+="-";
		
		stnum.clear();
		stnum << tS;	// current timestep is added to the name of the graphic outputfile
		stnum >> snum;
		
		procSS+= snum;	// procSS is a string process-numStep
		strcat(outGraphFile,".pvd");
		fGraph.open(outGraphFile);
		fGraph << "<?xml version=\"1.0\"?>\n";
		fGraph << "<VTKFile type=\"Collection\" version=\"0.1\" byte_order=\"LittleEndian\">\n";
		fGraph << "	<Collection>\n";
		fGraph.close();
		strcat(dirName,"-parts");
        mkdir(dirName,S_IRWXU);
	}
	if(!gFormat.compare("oogl"))	strcat(outGraphFile,".oog");
}

}

// ***

void Process::closeGraphOutFile(const string& gFormat, const int& tS) 
{
	if(!gFormat.empty())
	{
		if(!gFormat.compare("vtu"))	
		{
			fGraph.open(outGraphFile,ios::app);
			fGraph << "	</Collection>\n";
			fGraph << "</VTKFile>\n";
			fGraph.close();
		}
	}
}

// ***

void Process::closeFileDofs(void) 
{
	if(fDofs.is_open())		fDofs.close();
	if(fVelDofs.is_open())	fVelDofs.close();
}

// ***

void Process::closeFileStress(void) 
{
	if(fStress.is_open())		fStress.close();
}

// ***

void Process::setGravity(double g, int gd) 
{
	grav = g;
	gravDir = gd;
}

// ***

void Process::calcEner(void)
{
	Ekin = 0.0;
	Epot = 0.0;
	Econ = 0.0;
//      Eext = 0.0;
//      Dint = 0.0;
//	Eamb = 0.0;
	
	rootSys->calcEner(Ekin,Epot,Econ,Eext,Eamb,Dint,grav,gravDir,q,dq);
}

// ***

void Process::calcEntropy(void)
{
	Entropy = 0.0;
//	EntropyAmb = 0.0;
	
	rootSys->calcEntropy(Entropy, EntropyAmb);
}

// ***

void Process::calcAM(void)
{
	setZero(AM);
	rootSys->calcAM(AM,q,dq);
}

// ***

void Process::calcCon(double& time)
{
	// vCon has three elements: constraint and its two time derivatives
	setZero(vCon);
	rootSys->calcCon(vCon,q,dq,ddq,time);
}

// ***

void Process::sendIniPosToPoint(string& sName, string& iniPointName, vector<double>& iniValues)
{
	// Check if the point name has changed, in order to send it to the right system
	if (mMecSystems[sName]->checkIfPointChanged(iniPointName)==true)
	{
		iniPointName = mMecSystems[sName]->mChangedPoints[iniPointName];
		sName = getStringHead(iniPointName);
	}
	
	getPTSystem(sName)->setIniPosPoint(iniPointName,iniValues);
}

// ***

void Process::sendIniPosToPointSet(string& sName, string& iniPointSetName, vector<double>& iniValues)
{
	mMecSystems[sName]->setIniPosPointSet(iniPointSetName,iniValues);
}

// ***

void Process::sendIniVelToPoint(string& sName, string& iniPointName, vector<double>& iniValues)
{
	// Check if the point name has changed, in order to send it to the right system
	if (mMecSystems[sName]->checkIfPointChanged(iniPointName)==true)
	{
		iniPointName = mMecSystems[sName]->mChangedPoints[iniPointName];
		sName = getStringHead(iniPointName);
	}
//			cout << "sending initial velocity to point ... " << iniPointName << endl;
	getPTSystem(sName)->setIniVelPoint(iniPointName,iniValues);
//	mMecSystems[sName]->setIniVelPoint(iniPointName,iniValues);
}

// ***

void Process::sendIniVelToPointSet(string& sName, string& iniPointSetName, vector<double>& iniValues)
{
	mMecSystems[sName]->setIniVelPointSet(iniPointSetName,iniValues);
}

// ***

void Process::sendIniVelToBody(string& sName, string& bName, string& iniPointName, vector<double>& v0, vector<double>& Omega)
{
	mMecSystems[sName]->setIniVelBody(bName,iniPointName,v0,Omega);
}

// ***

void Process::setInitialConditions(double& h)
{
	if(rootSys!=0)	rootSys->setInitialConditions(q,dq);	/** Recursive from root system */
	
	// Imposed displacements impose initial velocities (for dynamics)
	//double tMinus = 0.0;
	int i;
	for (i=0; i<vImposedDofs.size(); ++i)
	{
		if (h!=0.0)
		{
			//dq[vImposedDofs[i]] =  ( mHistories[vImposedDofsHis[i]]->value(h) - mHistories[vImposedDofsHis[i]]->value(tMinus) )/ h;		
			dq[vImposedDofs[i]] =  ( mHistories[vImposedDofsHis[i]]->value(h) )/ h;		
			//cout << "dq[] = " << dq[vImposedDofs[i]] << endl;
		}
	}	
}

// ***

void Process::update(const int& inType)
{
	if(rootSys!=0)	rootSys->update(inType, qConv, dqConv);	/** Recursive from root system */
}

// ***

void Process::setInertiaDofNumbers()
{
	// Recursive from root system

	numDofsInertia = 0;	// Accumulates the dofs with inertia
	
	// WARNING: Inertia dofs are only at points at this moment, not at components
	
	// Set global numbers for translational and rotational dofs
	// numNoInertia collects the global dof numbers without inertia
	if(rootSys!=0)	numDofsInertia = rootSys->setInertiaNumbers(0,numNoInertia);			
	
	if (num==rankk) cout << "	--> Process " << num << " has " << numDofsInertia << " dofs with inertia \n";
//	if (num==rankk) cout << "	--> Process " << num << " has " << numNoInertia.size() << " dofs without inertia\n";
//	if (num==rankk) cout << "	numNoInertia es" << endl;
//	printSTLVector(numNoInertia);
}

// ***

void Process::createConvergedVectors()
{	
	if (qConv.size() == 0){qConv.resize(numDofs);}
	if (dqConv.size() == 0){dqConv.resize(numDofs);}
}

// ***

void Process::setQConv(void)
{
//	cout << "dyflag es " << dynFlag << endl;
	qConv = q;
	dqConv = dq;	
//	cout << "qConv es " << endl;
//	printSTLVector(qConv);
}

// ***

void Process::removeFromVector(vector<double>& vBig, vector<int>& vPos, vector<double>& vSmall)
{
	// Fills vSmall from vBig, based on the entries that are NOT specified in vPos
	// vPos must be ordered
	if (vPos.size()!=0)
	{
		int i;
		int cont = 0;
		for(i=0; i<vBig.size(); ++i)
		{
			if( vPos[cont]!=i)
			{
				vSmall.push_back( vBig[i] );
			}
			else
			{
				//			cout << "Me he saltado " << i << endl;
				cont+=1;
			}
		}
	}
	else vSmall = vBig;
}

// ***

void Process::restoreFromVector(vector<double>& vBig, vector<int>& vPos, vector<double>& vSmall)
{
	// Fills vBig from vSmall.  The entries specified in vPos are 0
	// vPos must be ordered
	
	if (vPos.size()!=0)
	{
		int i;
		int cont = 0;
		int cont2 = 0;
		for(i=0; i<vBig.size(); ++i)
		{
			if( vPos[cont]!=i)
			{
				vBig[i] = vSmall[cont2];
				cont2+= 1;
			}
			else 
			{
				vBig[i] = 0.0;
				cont+=1;
			}
		}
	}
	else vBig = vSmall;
}	

// ***

void Process::setTemperature(void) 
{
	rootSys->setTemperature(thetab);	/** Recursive from root system */
}

// ***

void Process::createTempOutFile()
{	
	stringstream stnum;
	string snum; 
	
	// Create if the output file does not exist
//	if ( !strcmp(outTempFile,"") )
	{
		stnum << num;
		stnum >> snum;
		
		strcpy(outTempFile,"Temp");
		strcat(outTempFile,snum.c_str());
		
		stnum.clear();
		
		stnum << rankk;
		stnum >> snum;
		
		strcat(outTempFile,".dat");
		
		fTemp.open(outTempFile);
		fTemp.close();
		fTemp.open(outTempFile,ios::app);
	}
}

// ***

void Process::fillVTempOut(const string& namein)
{
	fTemp << "#" << namein << endl;
	vTempOut.push_back(namein);	
}

// ***

void Process::writeFileTemp(const double& time)
{
	fTemp << time;
/*
	for (i=0; i<vDofsOut.size(); ++i)
	{
		fDofs << "\t" << q[ vDofsOut[i] ];
	}
	fDofs << "\n";
*/	
	string sN;
	for (int i=0; i<vTempOut.size(); ++i)
	{
		sN = getStringHead(vTempOut[i]);
		//cout << " a sistema " << sN << endl;
		getPTSystem(sN)->writeFileTemp(time,fTemp,vTempOut[i]);
	}
//	rootSys->writeFileTemp(time,fTemp);	/** Recursive from root system */
	fTemp << endl;
}

// ***

void Process::closeFileTemp(void) 
{
	if(fTemp.is_open())		fTemp.close();
}

// ***

void Process::sendIniPosToBody(string& sName, string& iniJointName, vector<double>& iniValues)
{
	mMecSystems[sName]->setIniPosBody(iniJointName,iniValues);
}

// ***

void Process::createIterOutFile()
{	
	stringstream stnum;
	string snum; 
	
	// Create if the output file does not exist
	//	if ( !strcmp(outTempFile,"") )
	stnum << num;
	stnum >> snum;
		
	strcpy(outIterFile,"Iter");
	strcat(outIterFile,snum.c_str());
	
	stnum.clear();
		
	stnum << rankk;
	stnum >> snum;
		
	strcat(outIterFile,".dat");
		
	fIter.open(outIterFile);
	fIter.close();
	fIter.open(outIterFile,ios::app);
}

// ***

void Process::writeFileIter(const double& time, const int& it)
{
	fIter << time << " " << it << endl;
}

// ***

void Process::closeFileIter(void) 
{
	if(fIter.is_open())		fIter.close();
}

// ***

void Process::setQatComponents(void) 
{
	rootSys->setQatComponents(q,dq);
}

// ***

void Process::updateLagMul(void) 
{
	rootSys->updateLagMul();
}

// ***

void Process::checkConvLagMul(bool& c) 
{
	rootSys->checkConvLagMul(c);
}

// ***

void Process::projectVel(double& projPenal, int facMode) 
{
	int i;
	int j;
	
	// Projection is performed as
	// (A + alpha*Phiq*Phiq)*qdot = A*qdot_ori
	
	M->copyMatrix(pMatrix);	
/*	
	cout << "M in projectVel" << endl;
	for (i=0; i<M->rows(); ++i)
	{
		for (j=0; j<M->rows(); ++j)
		{
			cout << M->readElement(i, j) << " ";
		}
		cout << endl;
	}
	cout << endl;
*/	
		
	vector<double> dqAux;
	
	removeFromVector(dq, numNoInertia, dqAux);
	
	// projVel = A*qdot_ori
	for (i=0; i<pMatrix->rows(); ++i)
	{
		projVel[i] = 0.0;
		for(j=0; j<pMatrix->rows(); ++j)
		{
			projVel[i]+= pMatrix->readElement(i,j)*dqAux[j];
		}
	}
	/*
	cout << "Original projection matrix A is " << endl;
	for (i=0; i<pMatrix->rows(); ++i)
	{
		for (j=0; j<pMatrix->rows(); ++j)
		{
			cout << pMatrix->readElement(i, j) << " ";
		}
		cout << endl;
	}
	cout << endl;
	 */
	
	// A + alpha*Phiq*Phiq
	rootSys->calcProjMatrix(*pMatrix, projPenal);

	/*
	cout << "Final projection matrix (A + penalty*Phiq*Phiq) is " << endl;
	for (i=0; i<pMatrix->rows(); ++i)
	{
		for (j=0; j<pMatrix->rows(); ++j)
		{
			cout << pMatrix->readElement(i, j) << " ";
		}
		cout << endl;
	}
	cout << endl;
	 */
	
	// Solve projection
	/*
	cout << "projVel antes is " << endl;
	printSTLVector(projVel);	
	*/
	pMatrix->factorize(facMode);	
	pMatrix->subsSolve(projVel);
	
	/*
	cout << "projVel despues is " << endl;
	printSTLVector(projVel);
	
	
	cout << "dq before... " << endl;
	printSTLVector(dq);
	 */
	
	// Substitution in the velocity vector	
	restoreFromVector(dq, numNoInertia, projVel);
	
	/*
	cout << "dq after.... " << endl;
	//cout << lTwoNorm(dq) << endl;
	printSTLVector(dq);
	 */
}

// ***

void Process::setPatternInertia(sparseMatrix* A)
{	
	A->sparsePattern(iaInertia,jaInertia);
	A->initMatrices();
}

// ***

double Process::calcProjPenalty(double& deltaProjTol, string& autoProjType)
{
	int i;
	int j;
	int k;
	vector<double> dotPhiAst;
	double delpRes;
//	double deltaExtreme = 0.0;
	double penal;
	double deltaEinf = 0.0;
	double dotPhiAstNorm = 0.0;
	
	// r is the number of constraint equations
	
	// Fill Phiq
	// Phiq is the constraint gradient.  r x numDofs
	Phiq.clear();
	rootSys->fillPhiq(Phiq,numDofs);
	
	//cout << "Phiq.size is " << Phiq.size() << endl;
	//cout << "Phiq[0].size is " << Phiq[0].size() << endl;
//	cout << "Phiq is " << endl;
//	printSTLMatrix(Phiq);
	//cout << "norma de Phiq[0] es " << lTwoNorm(Phiq[0]) << endl;
//	cout << "norma de Phiq[5] es " << lTwoNorm(Phiq[5]) << endl;
//	cout << "norma de Phiq[11] es " << lTwoNorm(Phiq[11]) << endl;
//	cout << "Phiq[5] es " << endl;
//	printSTLVector(Phiq[5]);
	
	auxGamma = Phiq;
	
	// Gamma is r x r
	if (Gamma.size()==0)
	{
		if (Phiq.size()!=0)
		{
			Gamma.resize(Phiq.size());
			for (i=0; i<Phiq.size(); ++i)	Gamma[i].resize(Phiq.size());
		}
	}
	// M*x = PhiqT
	// x es (r x numDofs)
	
	MfacProj->subsSolveTrans(auxGamma);
	
	// Warning: auxGamma is transposed (r x numDofs)
	
//	cout << "La solucion es " << endl;
//	printSTLMatrix(auxGamma);
	
	// Gamma = Phiq*x
	for (i=0; i<Phiq.size(); ++i)
	{
		for(j=0; j<Phiq.size(); ++j)
		{
			Gamma[i][j] = 0.0;
			for(k=0; k<Phiq[0].size(); ++k)
			{
				
				Gamma[i][j] += Phiq[i][k]*auxGamma[j][k];
			}
		}
	}
//	cout << "Gamma is " << endl;
//	printSTLMatrix(Gamma);
	
	GammaFac = Gamma;
	
	// dotPhiAst is the velocity constraint vector before projection
	// dotPhiAst = Phiq*dq
	dotPhiAst.resize(Phiq.size());
	for (i=0; i<Phiq.size(); ++i)
	{
		dotPhiAst[i] = 0.0;
		for(j=0; j<numDofs; ++j)
		{
			dotPhiAst[i]+= Phiq[i][j]*dq[j];
		}
	}
	
//	cout << "dotPhiAst is " << endl;
//	printSTLVector(dotPhiAst);

	dotPhiAstNorm = lTwoNorm(dotPhiAst);
	
//	cout << "AutProjFlag = " << autoProjFlag << endl;
//	cout << "AutProjType = " << autoProjType << endl;
//	cout << "dotPhiAstNorm is " << dotPhiAstNorm <<  endl;
	if (autoProjFlag==1)	// Energy criterions
	{
		// Maximum energy dissipated in projection
		// delEpInf = (1/2) dotPhiAstT * Gamma-1 * dotPhiAst
		vector<double> z;
		z = dotPhiAst;
		lu(GammaFac);
		forSub(GammaFac,z);
		backSub(GammaFac,z);
		for (i=0; i<z.size(); ++i)
		{
			deltaEinf += 0.5*dotPhiAst[i]*z[i];
		}
//		cout << "deltaEinf = " << deltaEinf << endl;

		// Below the tolerance, use a large penalty value
	}
	
	
	if ( autoProjFlag==1 && deltaProjTol==0.0)			penal = 0.0;
	else if ( autoProjFlag==1 && deltaProjTol>0.0 && deltaEinf<=deltaProjTol)	penal = projPenaltyMax;	// Energy already below tolerance. Maximum alpha
	else if ( autoProjFlag==2 && dotPhiAstNorm<=deltaProjTol)		penal = 0.0; 
	else if (!autoProjType.compare("Exact") ) // Bisection method for 'exactly' solving DeltaEp or DotPhi  													// alpha must be calculated.
	{
		// Bisection between alpha=0 and alpha=1.0e12
		double al1 = 0.0;
		double al2 = projPenaltyMax;
		double fal1;
		double fal2;
		double atol = 1.0e-6;
		double mida;
		double fmida;
		
		int maxiter = ceil(log2 (al2-al1) - log2 (2*atol));
//		cout << "maxiter = " << maxiter << endl;
		
		//cout << "al1 = " << al1 << endl;
		//cout << "al2 = " << al2 << endl;
		if (autoProjFlag==1)	// Energy criterion
		{
			fal1 = calcDeltaEpRes(al1, Gamma, dotPhiAst, deltaProjTol);
			fal2 = calcDeltaEpRes(al2, Gamma, dotPhiAst, deltaProjTol);
		}
		else if (autoProjFlag==2)	// Velocity constraint criterion
		{
			fal1 = calcDotPhiRes(al1, Gamma, dotPhiAst, deltaProjTol);
			fal2 = calcDotPhiRes(al2, Gamma, dotPhiAst, deltaProjTol);
		}
		//cout << "fal1 = " << fal1 << endl;
		//cout << "fal2 = " << fal2 << endl;
	
		for (i=0; i<maxiter; ++i)
		{
			mida = (al1+al2)/2.0;
			if (autoProjFlag==1)		fmida = calcDeltaEpRes(mida, Gamma, dotPhiAst, deltaProjTol);
			else if (autoProjFlag==2)	fmida = calcDotPhiRes(mida, Gamma, dotPhiAst, deltaProjTol);
//			cout << "mida = " << mida << endl;
//			cout << "fmida = " << fmida << endl;
			if (fal1*fmida < 0.0)
			{
				al2 = mida;
				fal2 = fmida;
			}
			else 
			{
				al1 = mida;
				fal1 = fmida;
			}
		}
		penal = (al1+al2)/2.0;
		//cout << "******* Este es el fmid " << calcDotPhiRes(mida, Gamma, dotPhiAst, deltaProjTol) << endl;
		//cout << "******* Este es el penal " << penal << endl;
//		cout << "Gamma is " << endl; 
//		printSTLMatrix(Gamma);
//				cout << "******* estaba en ener exact " << endl;
	}
	else if (autoProjFlag==1 && !autoProjType.compare("Exp") ) // Exponential function approximation for DeltaEp
	{
	  double r, s;
	  
	  r = -deltaEinf;
	  
	  s = dotPhiAstNorm*dotPhiAstNorm/r;
		  
	  penal = 1.0/s * log(1.0 + deltaProjTol/r);
	}
	else if (autoProjFlag==1 && !autoProjType.compare("Rat") ) // Rational function approximation for DeltaEp
	{ 
	   penal = deltaProjTol/(dotPhiAstNorm*dotPhiAstNorm*(1 - deltaProjTol/deltaEinf));
	}
	else if (autoProjFlag==1 && !autoProjType.compare("Quad") ) // Quadratic function approximation for DeltaEp
	{ 
		double a = 0.0;
		double b = dotPhiAstNorm*dotPhiAstNorm;
		double sol1;
		double sol2;
		for (i=0; i<dotPhiAst.size(); ++i)
		{
			for(j=0; j<dotPhiAst.size(); ++j)
			{
				a -= dotPhiAst[i]*Gamma[i][j]*dotPhiAst[j];
			}
		}
		a *= 3.0;
		/*
		cout << "a = " << a << endl;
		cout << "b = " << b << endl;
		cout << "deltaProjTol = " << deltaProjTol << endl;
		cout << "inside sqrt = " << b*b+4.0*a*deltaProjTol << endl;
		cout << "DeltaEinf = " << deltaEinf << endl;
		*/
		sol1 = (-b + sqrt(b*b+4.0*a*deltaProjTol))/2.0/a;
		sol2 = (-b - sqrt(b*b+4.0*a*deltaProjTol))/2.0/a;
		//cout << "sol1 = " << sol1 << endl;
		//cout << "sol2 = " << sol2 << endl;
		if (sol1<=sol2)	penal = sol1; else penal = sol2;
	}
	else if (autoProjFlag==1 && !autoProjType.compare("Linear") ) // Linear function approximation for DeltaEp
	{
		penal = deltaProjTol/dotPhiAstNorm/dotPhiAstNorm;
	}
	else if (autoProjFlag==2 && !autoProjType.compare("Exp") ) // Exponential function approximation for DotPhi
	{
	  double f = 0.0;
	  
	  for (i=0; i<dotPhiAst.size(); ++i)
	  {
	    for(j=0; j<dotPhiAst.size(); ++j)
	    {
		    f -= dotPhiAst[i]*Gamma[i][j]*dotPhiAst[j];
	    }
	  }
	  f = f/dotPhiAstNorm/dotPhiAstNorm;
	  penal = 1.0/f * log(deltaProjTol/dotPhiAstNorm);
	}
	else if (autoProjFlag==2 && !autoProjType.compare("Rat") ) // Rational function approximation for Dotphi 
	{
		double a = 0.0;
	  
		for (i=0; i<dotPhiAst.size(); ++i)
		{
			for(j=0; j<dotPhiAst.size(); ++j)
			{
				a += (dotPhiAst[i]*Gamma[i][j]*dotPhiAst[j]);
			}
		}
		a = dotPhiAstNorm*dotPhiAstNorm/a;	
		penal = sqrt(dotPhiAst.size())*a*dotPhiAstNorm/deltaProjTol - a;
	}
	else if (autoProjFlag==2 && !autoProjType.compare("Inv"))
    {
        vector<double> g;
        g = dotPhiAst;
        lu(GammaFac);
        forSub(GammaFac,g);
        backSub(GammaFac,g);
        
        penal = lTwoNorm(g)/ deltaProjTol;
    }	
	return penal;
}

// ***

double Process::calcDeltaEpRes(double alpha, vector<vector<double> >& Gam, vector<double>& dP, double deltaEpTol)
{
	int i;
	int j;

	vector<vector<double> > onePlusAlGam;
	vector<double> xPlusdP;
	double aux = 0.0;
	
	
	onePlusAlGam = Gam;
	xPlusdP = dP;

	// onePlusAlGam = 1 + alpha*Gam
	for (i=0; i<onePlusAlGam.size(); ++i)
	{
		for(j=0; j<onePlusAlGam.size(); ++j)
		{
			onePlusAlGam[i][j]*= alpha;
		}
		onePlusAlGam[i][i] += 1.0;
	}

	lu(onePlusAlGam);
	
	// onePlusAlGam * x = dP
	forSub(onePlusAlGam,xPlusdP);
	backSub(onePlusAlGam,xPlusdP);
	
//	cout << "x = " << endl;
//	printSTLVector(xPlusdP);
	
	for (i=0; i<xPlusdP.size(); ++i)
	{
		xPlusdP[i]+= dP[i];
	}

	// onePlusAlGam * y = (x + dP)
	forSub(onePlusAlGam,xPlusdP);
	backSub(onePlusAlGam,xPlusdP);
	
//	cout << "y = " << endl;
//	printSTLVector(xPlusdP);

	
	for (i=0; i<dP.size(); ++i)
	{
		aux+= dP[i]*xPlusdP[i];
	}
	aux*= alpha/2.0;
	aux-= deltaEpTol;
	
	return aux;
}

// ***

double Process::calcDotPhiRes(double alpha, vector<vector<double> >& Gam, vector<double>& dP, double deltaTol)
{
	int i;
	int j;
	
	vector<vector<double> > onePlusAlGam;
	vector<double> xPlusdP;
	
	onePlusAlGam = Gam;
	xPlusdP = dP;
	
	// onePlusAlGam = 1 + alpha*Gam
	for (i=0; i<onePlusAlGam.size(); ++i)
	{
		for(j=0; j<onePlusAlGam.size(); ++j)
		{
			onePlusAlGam[i][j]*= alpha;
		}
		onePlusAlGam[i][i] += 1.0;
	}
	
	lu(onePlusAlGam);
	
	// onePlusAlGam * x = dP
	forSub(onePlusAlGam,xPlusdP);
	backSub(onePlusAlGam,xPlusdP);
	
//	cout << "alpha***** = " << alpha << endl;
	//cout << "xPlus = " << endl;
	//printSTLVector(xPlusdP);

//	cout << "xPlusdP ***** = " << lTwoNorm(xPlusdP) << endl;
	return lTwoNorm(xPlusdP)-deltaTol;
}

// ***

void Process::setSurfaces(const map<string, Surface*, less<string> >& mSurf)
{
	mSurfaces = mSurf;
}

// ***

void Process::setSurfAtBodies()
{
	// Recursive from root system
	if(rootSys!=0)	rootSys->setSurfAtBodies();
}

// ***

void Process::createContacts()
{
	// Recursive from root system
	if(rootSys!=0)	rootSys->createContacts();
}

// ***

void Process::addNoDofsPoint(string& pname, string& stDofs, double& endtime)
{
	mPointsNoDofs.insert(pair<string,string>(pname,stDofs));
	vPointsEndTimes.push_back(endtime);
}

// ***

void Process::setNoDofs()
{
	int i;
	multimap<string, string, less<string> >::iterator it;
	string auxs;
	vector<int> vauxi;
	int cont = 0;
	
	//cout << "vPointsEndTimes = " << endl;
	//printSTLVector(vPointsEndTimes);
	
	for(it=mPointsNoDofs.begin(); it!=mPointsNoDofs.end(); it++)
	{
		auxs = it->second;
		vauxi = getPTPoint(it->first)->getGlobalNumbers();
		
		if (!auxs.compare("all"))
		{
			for (i=0; i<vauxi.size(); ++i)
			{
				vNoDofs.push_back(vauxi[i]);
				vNoDofsEndTimes.push_back(vPointsEndTimes[cont]);
			}
		}		
		else if (!auxs.compare("123"))
		{
			vNoDofs.push_back(vauxi[0]);
			vNoDofs.push_back(vauxi[1]);
			vNoDofs.push_back(vauxi[2]);
			vNoDofsEndTimes.push_back(vPointsEndTimes[cont]);
			vNoDofsEndTimes.push_back(vPointsEndTimes[cont]);
			vNoDofsEndTimes.push_back(vPointsEndTimes[cont]);
		}
		else if (!auxs.compare("12"))
		{
			vNoDofs.push_back(vauxi[0]);
			vNoDofs.push_back(vauxi[1]);
			vNoDofsEndTimes.push_back(vPointsEndTimes[cont]);
			vNoDofsEndTimes.push_back(vPointsEndTimes[cont]);
		}
		else if (!auxs.compare("23"))
		{
			vNoDofs.push_back(vauxi[1]);
			vNoDofs.push_back(vauxi[2]);
			vNoDofsEndTimes.push_back(vPointsEndTimes[cont]);
			vNoDofsEndTimes.push_back(vPointsEndTimes[cont]);
		}
		else if (!auxs.compare("13"))
		{
			vNoDofs.push_back(vauxi[0]);
			vNoDofs.push_back(vauxi[2]);
			vNoDofsEndTimes.push_back(vPointsEndTimes[cont]);
			vNoDofsEndTimes.push_back(vPointsEndTimes[cont]);
		}
		else if (!auxs.compare("1"))
		{
			vNoDofs.push_back(vauxi[0]);
			vNoDofsEndTimes.push_back(vPointsEndTimes[cont]);
		}
		else if (!auxs.compare("2"))
		{
			vNoDofs.push_back(vauxi[1]);
			vNoDofsEndTimes.push_back(vPointsEndTimes[cont]);
		}
		else if (!auxs.compare("3"))
		{
			vNoDofs.push_back(vauxi[2]);
			vNoDofsEndTimes.push_back(vPointsEndTimes[cont]);
		}	
		vauxi.clear();
		cont +=1;
	}
	/*
	cout << "vNoDofs is " << endl;
	printSTLVector(vNoDofs);
	cout << "vNoDofsEndTimes is " << endl;
	printSTLVector(vNoDofsEndTimes);
	 */
}

// ***

void Process::createLumpedMass()
{	
	vLumpedMass.resize(numDofsInertia);
	rootSys->fillMassMatrix(vLumpedMass,q,0.0);
/*	
	int i;
	for(i=0; i<vInterfaces.size(); i++)	
	{
		vInterfaces[i].setLumpedMass(vLumpedMass);
	}	
*/
}

// ***

void Process::synchronizeLumpedMass()
{	
	int i;
	for(i=0; i<vInterfaces.size(); i++)	
	{
		vInterfaces[i].synchronizeMass(vLumpedMass);
	}
}

// ***

void Process::repairChangedPointNames()
{
	rootSys->repairChangedPointNames();
}

// ***

void Process::fillMSmooth(double& time)
{
	if (matType==1)
	{
		MSmooth = new sparseMatrixSLU();
	}
	if (matType==2)
	{
		MSmooth = new sparseMatrixUMF();
	}
	if (ia.size()==0)			setSparsePattern();
	if (iaInertia.size()==0)	setSparsePatternInertia();

	setPatternInertia(MSmooth);

	rootSys->fillMSmooth(*MSmooth,q,time);
	
	//cout << "MSmooth es " << endl;
	//printSparseMatrix(MSmooth);
	
	int kk1 = 0;
	int kk2 = 0;
	
	MSmooth->factorize(0);
}

// ***

void Process::calcSmoothStress(double& time)
{
	// FSmooth has 2 rows
	// First: sigma_x, sigma_y, sigma_z
	// Second: sigma_xy, sigma_xz, sigma_yz
	int i;
	if (FSmooth.size()==0)	
	{
		//FSmooth.resize(4);
		FSmooth.resize(2);
		//for (i=0; i<4; ++i)	FSmooth[i].resize(numDofsInertia);
		for (i=0; i<2; ++i)	FSmooth[i].resize(numDofsInertia);
	}
	else setZero(FSmooth);
	
	if (FSmoothDerived.size()==0)	
	{
		//FSmooth.resize(4);
		FSmoothDerived.resize(2);
		//for (i=0; i<4; ++i)	FSmooth[i].resize(numDofsInertia);
		for (i=0; i<2; ++i)	FSmoothDerived[i].resize(numDofsInertia);
	}
	else setZero(FSmoothDerived);

	rootSys->fillFSmooth(FSmooth,q,dq,time);

	//cout << "FSSSSSmmm [0][0]" << FSmooth[0][0] <<  endl;
	//printSTLMatrix(FSmooth);
	
	MSmooth->subsSolveTrans(FSmooth);
	
	// Computation of principal and Von Mises stresses based on the smoothed stresses
	if (PSigma.size()==0)	PSigma.resize(dimension);
	if (SigmaMat.size()==0)	SigmaMat.resize(dimension,dimension);
	
	int j;
	int k;
	int nnum = (int) numDofsInertia/dimension;
	//cout << "nnum es " << nnum << endl;
	// Loop on the nodes with inertia
	for (k=0; k<nnum; ++k)
	{
		// Diagonal terms
		for (i=0; i<dimension; ++i)	SigmaMat(i,i) = FSmooth[0][dimension*k+i];
	
		// Off-diagonal
		if (dimension==2)
		{
			SigmaMat(0,1) = FSmooth[1][dimension*k];
			SigmaMat(1,0) = SigmaMat(0,1);
		}
		if (dimension==3)
		{
			SigmaMat(0,1) = FSmooth[1][dimension*k];
			SigmaMat(0,2) = FSmooth[1][dimension*k+1];
			SigmaMat(1,2) = FSmooth[1][dimension*k+2];
			SigmaMat(1,0) = SigmaMat(0,1);
			SigmaMat(2,0) = SigmaMat(0,2);
			SigmaMat(2,1) = SigmaMat(1,2);			
		}
		// Ordered eigenvalues (largest first)
		SelfAdjointEigenSolver<MatrixXd> eigensolver(SigmaMat);
		for (i=0; i<dimension; ++i)	PSigma[i] = eigensolver.eigenvalues()[i];		
		
		sort(PSigma.begin(), PSigma.end());
		
		reverse(PSigma.begin(), PSigma.end());
		
		//for (i=0; i<dimension; ++i)	FSmooth[2][dimension*k+i] = PSigma[i];
		for (i=0; i<dimension; ++i)	FSmoothDerived[0][dimension*k+i] = PSigma[i];
		
		// Von Mises
		//FSmooth[3][dimension*k] = sqrt( ( (PSigma[0]-PSigma[1])*(PSigma[0]-PSigma[1]) +  (PSigma[1]-PSigma[2])*(PSigma[1]-PSigma[2]) +  (PSigma[2]-PSigma[0])*(PSigma[2]-PSigma[0]) )/2.0 );
		FSmoothDerived[1][dimension*k] = sqrt( ( (PSigma[0]-PSigma[1])*(PSigma[0]-PSigma[1]) +  (PSigma[1]-PSigma[2])*(PSigma[1]-PSigma[2]) +  (PSigma[2]-PSigma[0])*(PSigma[2]-PSigma[0]) )/2.0 );
		// Other stress measurements, TO BE DONE
		/*
		FSmooth[3][dimension*k+1] = FSmooth[3][dimension*k];
		FSmooth[3][dimension*k+2] = FSmooth[3][dimension*k];		
		*/
		FSmoothDerived[1][dimension*k+1] = FSmoothDerived[1][dimension*k];
		FSmoothDerived[1][dimension*k+2] = FSmoothDerived[1][dimension*k];		
	}
}

// ***

void Process::addImposedDispPoint(string& pname ,string& dofnum, char* fn)
{
	vector<string> saux(2);
	string fileName = fn;
	saux[0] = dofnum;
	saux[1] = fileName;
	//cout << "Aqui he llegado......" << pname << " , " << dofnum << " , " << fn << endl;
	mPointsImposedDofs.insert(pair<string,vector<string> >(pname,saux));
}

// ***

void Process::addImposedDispPointSet(string& psetname ,string& dofnum, char* fn)
{
	int i;
	string sysName = getStringHead(psetname);
	vector<string> vName = mMecSystems[sysName]->pointSets[psetname];
	
	for(i=0; i<vName.size(); ++i)
	{
		addImposedDispPoint(vName[i], dofnum, fn);
	}	
}

// ***

void Process::setImposedDofs()
{
	int i;
	//multimap<string, string, less<string> >::iterator it;
	multimap<string, vector<string>, less<string> >::iterator it;
	string auxs0;
	string auxs1;
	vector<int> vauxi;
	
	for(it=mPointsImposedDofs.begin(); it!=mPointsImposedDofs.end(); it++)
	{
		//auxs = it->second;
		auxs0 = it->second[0];
		auxs1 = it->second[1];

		vauxi = getPTPoint(it->first)->getGlobalNumbers();
		
		if (!auxs0.compare("all"))
		{
			for (i=0; i<vauxi.size(); ++i)
			{
				vImposedDofs.push_back(vauxi[i]);
				vImposedDofsHis.push_back(auxs1);
				//mImposedDofs[vauxi[i]] = auxs1;
			}
		}		
		else if (!auxs0.compare("123"))
		{
			
			vImposedDofs.push_back(vauxi[0]);
			vImposedDofs.push_back(vauxi[1]);
			vImposedDofs.push_back(vauxi[2]);
			vImposedDofsHis.push_back(auxs1);
			vImposedDofsHis.push_back(auxs1);
			vImposedDofsHis.push_back(auxs1);			
			/*
			mImposedDofs[vauxi[0]] = auxs1;
			mImposedDofs[vauxi[1]] = auxs1;
			mImposedDofs[vauxi[2]] = auxs1;
			 */
		}
		else if (!auxs0.compare("12"))
		{
			vImposedDofs.push_back(vauxi[0]);
			vImposedDofs.push_back(vauxi[1]);
			vImposedDofsHis.push_back(auxs1);
			vImposedDofsHis.push_back(auxs1);
			/*
			mImposedDofs[vauxi[0]] = auxs1;
			mImposedDofs[vauxi[1]] = auxs1;
			*/
		}
		else if (!auxs0.compare("23"))
		{
			vImposedDofs.push_back(vauxi[1]);
			vImposedDofs.push_back(vauxi[2]);
			vImposedDofsHis.push_back(auxs1);
			vImposedDofsHis.push_back(auxs1);
			/*
			mImposedDofs[vauxi[1]] = auxs1;
			mImposedDofs[vauxi[2]] = auxs1;
			*/
		}
		else if (!auxs0.compare("13"))
		{
			
			vImposedDofs.push_back(vauxi[0]);
			vImposedDofs.push_back(vauxi[2]);
			vImposedDofsHis.push_back(auxs1);
			vImposedDofsHis.push_back(auxs1);
			/*
			mImposedDofs[vauxi[0]] = auxs1;
			mImposedDofs[vauxi[2]] = auxs1;
			*/
		}
		else if (!auxs0.compare("1"))
		{
			vImposedDofs.push_back(vauxi[0]);
			vImposedDofsHis.push_back(auxs1);
			//mImposedDofs[vauxi[0]] = auxs1;
		}
		else if (!auxs0.compare("2"))
		{
			vImposedDofs.push_back(vauxi[1]);
			vImposedDofsHis.push_back(auxs1);
			//mImposedDofs[vauxi[1]] = auxs1;
		}
		else if (!auxs0.compare("3"))
		{
			vImposedDofs.push_back(vauxi[2]);
			vImposedDofsHis.push_back(auxs1);
			//mImposedDofs[vauxi[2]] = auxs1;
		}	
		vauxi.clear();
	}
	/*
	map<int, string, less<int> >::iterator itkk;
	for(itkk=mImposedDofs.begin(); itkk!=mImposedDofs.end(); itkk++)
	{
		cout << "mImposedDofs[" << itkk->first << "] = " << itkk->second << endl;
	}
	 */
	/*
	cout << "vImposedDofs is " << endl;
	printSTLVector(vImposedDofs);	
	cout << "vImposedDofsHis is " << endl;
	printSTLVector(vImposedDofsHis);	
	*/
}

// ***

void Process::setHistories(map<string,History*,less<string> >& mHis)
{
	mHistories = mHis;
}

// ***

void Process::setQIni()
{
	qIni = q;
}

int Process::checkComponents()
{
	int state = 0;
	if(rootSys!=0)	state = rootSys->checkComponents(state);	/** Recursive from root system */
	//cout << "process:state is " << state << endl;
	return state;
}

